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1. AN OVERVIEW OF THE ROTATING SPACECRAFT-TANK 
CONVECTION PROBLEM 

E. Dale Martin 




SUMMARY 


This introductory article provides the background for the succeeding 
articles comprising a study of convection in the tanks of a rotating space- 
craft. The discussion relates the analysis of mixing of a stratified fluid in 
a supercritical cryogenic state (e.g. , oxygen, in which mixing fans may not be 
used) to spacecraft such as Apollo and orbiting space stations. 


INTRODUCTION 


This report describes a general study of convection and mixing of a strat- 
ified fluid in a rotating container and its application to the special problem 
of fluid heating and convection in a spacecraft tank. In a rotating container, 
spatial variations of temperature and density lead to "natural convection" 
under the action of the effective body forces induced by the rotation. Simul- 
taneously, time-dependent rotation of the container may induce "forced convec- 
tion" coupled with the effects of stratification (see ch. 2). 

The study reported here began with the need to analyze the oxygen storage 
system used in an Apollo spacecraft, as outlined in the next section. There- 
fore , although the study applies to spacecraft fluid storage systems in gen- 
eral (including, e.g., orbiting space stations), the following conditions are 
imposed (at certain stages in the study) in the specific application to an 
Apollo oxygen tank: the fluid considered is oxygen in a supercritical cryo- 

genic state, near the critical point (where the fluid is not distinguishable 
as being either a gas or a liquid) ; the system is free of gravitational 
forces ; and the container is rotating about a noncentral axis . 

In a near-critical thermodynamic state , a fluid such as oxygen may 
undergo a phenomenon known as "pressure collapse," which is a local rapid 
decrease in pressure resulting from the sudden mixing of a thermally strati- 
fied fluid. The potential pressure drop (also called "potential pressure 
decay" and "potential for pressure collapse") is defined as the difference 
between the pressure in a thermally stratified environment and the pressure 
that would be present if the stratified fluid were suddenly mixed adiabati- 
cally into a uniform state. Thus, because of the nature of the state rela- 
tions in the near-critical region (where the state is highly sensitive to 
changes in temperature), sudden mixing of a significantly stratified super- 
critical fluid results in a substantial pressure drop. From many state 
points , such an occurrence would drive the fluid into the undesirable two- 
phase state, where vapor and liquid can coexist. The sudden decrease in pres- 
sure may render inoperative the systems (such as fuel cells and life support) 
intended to be supplied by the oxygen-storage system. This phenomenon of 
pressure collapse can be prevented by limiting the stratification, which is 
done by keeping the fluid sufficiently mixed. (Slow mixing limits the 
possibility of sudden mixing. ) 

From the above description, one can surmise that results of calculations 
of flow of a near-critical fluid may be highly dependent on accurate 
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representation of the thermodynamic properties (cf. ch. k and 6). Near- 
critical fluids also have low viscosity coefficients and high thermal expan- 
sion coefficients (ref. 1, p. 660 ff) so there will he strong coupling between 
the natural convection and the forced convection due to time-dependent 
rotation (cf. ch. 2). 

The goals of the study reported here are to simulate - by numerical anal- 
ysis and computation - the fluid motion, the energv transfer, and the mixing 
of a fluid under the above conditions , and to determine in particular the 
effectiveness of the mixing resulting from rotation maneuvers in reducing the 
potential pressure drop in the cryogenic oxygen-storage system used in space- 
craft such as the Apollo. 


BACKGROUND FROM APOLLO 


In spacecraft such as Apollo, supercritical cryogenic storage of oxygen 
and hydrogen is highly desirable because a large mass of the fluid can be 
stored at very high density, and therefore in a small volume, at pressures low 
enough for reasonably lightweight tanks. 

The Apollo oxygen tanks are located in the service module as shown in 
figure 1.1 for Apollo 13. The tanks are spheres with inside diameters of 25 
in. and with distances from the spacecraft axis to the sphere centers of 3 ft 
(tank l) and about 5 ft (tank 2). With the assumption that the center of mass 
of the vehicle is on the central axis , these distances are then the "rotation 
arms" of the spheres when the vehicle rotates in space. Each tank is loaded 
initially with 330 lb of liquid oxygen. The oxygen is then heated by an 
internal electrical heater (see fig. 1.2) to completely vaporize the oxygen. 
Subsequently, the heater is operated periodically to maintain a system design 
pressure of 900 *35 psia. The maintenance of system pressure by the heater 
operation is required because pressure is changed both by "heat leak" (heat 
conduction into the tank) at the imperfectly insulated wall and by withdrawal 
of fluid from the tank. Pressure sensors switch on the heater when the pres- 
sure drops to the desired lower limit and switch off the heater at the upper 
pressure limit. These "heater cycles" and corresponding "pressure cycles" are 
discussed in chapter 6. 

If there were no convection of the fluid (in zero gravity) during heater 
operation, strong temperature and density gradients would develop because of 
the inefficiency of pure heat conduction for energy transfer. Each heater 
cycle results in increasingly severe gradients; therefore, a significant 
potential pressure drop would develop so that a small subsequent acceleration 
could result in pressure collapse. Adequate slow mixing of the fluid there- 
fore is required to prevent severe stratification. 

In manned space missions prior to Apollo l4, electrically driven fans 
inside the oxygen tanks (see fig. 1.2) were used to ensure sufficient circula- 
tion of the fluid. A series of incidents resulted in electric arcing of the 
fan wiring in the no. 2 oxygen tank of Apollo 13, with combustion of Teflon 
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wire insulation leading to failure of 
the conduit-entry plug* a sudden 
release of pressure from the tank* and 
ultimate abortion of the mission. 


Closeout cap 


Overboard 


Relief valve 


Figure 1.2.- Oxygen tank no. 2 internal components 
(Apollo 13). 


In June 1970* the Apollo 13 
Review Board recommended removal of 
the wiring and fan motors* and also 
suggested further investigation of the 
need for stirring in the cryogenic 
oxygen-storage system. In the 
redesign effort hy MSA Manned Space- 
craft Center* an Apollo Cryogenic 
Oxygen Tank Analysis Team was estab- 
lished, with one function being to 
analyze the stratification problem. 
Various investigators contributed to 
that team effort. With the stratifi- 
cation established as possibly crucial 
by the team, one problem was to deter- 
mine the adequacy of vehicle rotation 
and maneuvers to produce sufficient 
convection for the required mixing in 
The contributions by Ames Research Center 


the oxygen tanks , without fans . 
were based on the studies described in this report; preliminary results were 
reported at a meeting of the analysis team by B. Baldwin and Y. Sheaf fer in 
January 1971 and at the NASA-MSC Cryogenics Symposium by B. Baldwin* E. D. 
Martin, W. A. Reinhardt, and Y. Sheaffer in May 1971. 


APPROACH AND SCOPE OF THE STUDY 


The format 
as chapters * is 
overall subject 
fields . 


of this report* a collection of separate but related articles 
used because each article makes contributions both to the 
of the report and to general knowledge in separate specialized 


In chapter 2 , Martin and Baldwin develop a set of approximate equations 
for the Navier-Stokes description of fluid convection in a rotating system, 
including the influence both of effective buoyancy body forces due to tempera- 
ture and density stratifications and of arbitrary time-dependent rotation and 
acceleration of the tank. The equations represent a small-density- variation 
approximation that is valid under conditions of low relative velocity but with 
significant density and temperature variations. All relevant terms represent- 
ing effects of rotation and changes in rotation are included. For example, 
the analysis includes the Coriolis terms known to be significant in the three- 
dimensional flow, as shown by experiments with dye injected into a water- 
filled rotating sphere (unpublished results by J. F. Lands, Jr., and R. C. 
Ried, Jr., reported at a meeting of the Apollo Tank Analysis Team in January 
1971). The equations are valid for combined forced and contained natural con- 
vection, whereas previous treatments of convection in rotating systems had 
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dealt with either pure natural convection or pure forced convection. For sub- 
sequent computations, the two-dimensional square-tank problem is formulated. 

In chapter 3, Lomax and Bailey describe a numerical finite-differencing 
scheme for computing the convection of vorticity and energy in a two- 
dimensional rotating tank, starting with the partial differential equations 
derived in chapter 2. The highly efficient computational method employs at 
each time step a Buneman Poisson solver for the Poisson equation relating the 
stream function to vorticity, and uses MacCormack's method for the transport 
equations for vorticity and energy. 

Reinhardt then develops efficient and versatile methods for the numerical 
evaluation of thermodynamic properties of cryogenic oxygen from semiempirical 
equations (ch. k ) . For reasons discussed above, accurate thermodynamic- 
property representations are needed for computations in the near-critical 
region. In the efficient methods developed, an arbitrary choice of indepen- 
dent state variables is allowed. A procedure for the rapid evaluation of 
volume integrals of spatially dependent quantities is introduced that depends 
on temperature distribution functions. By this method, insignificant losses 
in accuracy are accepted in exchange for significant savings in computation 
time. 


In chapter 5, Baldwin, Reinhardt, and Sheaffer discuss how the results 
based on the small-density-variation approximation (with density considered 
only as a function of temperature and with pressure as a slowly varying 
parameter) can be interpreted to determine the thermodynamic quantities. The 
latter include the slowly varying average pressure, the actual density dis- 
tribution, and the potential pressure drop. The problem of evaluating the 
time-dependent thermodynamic state of stratified fluid is considered based on 
(l) van der Waals equations of state, and (2) more exact thermodynamic 
relations . 

In chapter 6, Baldwin, Reinhardt, and Sheaffer present results of the 
numerical simulation for studying the effectiveness of rotation reversal and 
spinup from rest on mixing stratified oxygen in the Apollo storage tanks. The 
computations are two-dimensional simulations based on the complete equations 
developed in chapter 2, with use of the methods from chapters k and 5 for the 
accurate thermodynamics and efficient evaluations , and with use of the compu- 
tational method described in chapter 3. The significant effects of the rota- 
tion maneuvers on the potential pressure drop are discussed. 

Each chapter concludes with a list of references. 


REFERENCE 


1. Ostrach, S.: Laminar Flows with Body Forces. Section F in Theory of 

Laminar Flows, F. K. Moore, ed. , Vol. IV of High Speed Aerodynamics and 
Jet Propulsion, Princeton Univ. Press, 1964, pp. 528-718. 
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2. FLOW EQUATIONS GOVERNING COMBINED FORCED AND NATURAL CONVECTION 
FOR APPLICATION IN A ROTATING TANK 


E. Dale Martin and Barrett S. Baldwin 



SUMMARY 


A set of simplified equations is derived that represents a small-density- 
variation approximation of the Navier-Stokes equations for combined forced and 
contained natural convection of a stratified fluid. These equations reduce to 
a generalization of the Boussinesq approximation in the special case of steady 
rotation with (variable) buoyancy body forces, but also properly represent the 
forced convection due to time-dependent rotation. For use in subsequent 
numerical computations, the problem of convection in a rotating square tank 
with temperature and density stratifications is formulated. 


INTRODUCTION 


In this chapter, equations for fluid convection in a tank with arbitrary 
time-dependent rotation are developed from the Navier-Stokes equations. Sim- 
plified forms of these equations are obtained for use in numerical computation 
of the fluid convection in noninertial coordinate systems. Conditions assumed 
to exist in the rotating tank are small temperature and density variations and 
small apparent body forces . 

To clarify the relationship of the present problem to other convection 
problems, consider the terminology of, and conditions for, various classes of 
convection. Two special cases of convection (flow) are forced convection and 
natural convection (see Prandtl, ref. 1, pp. 396 and 4l2). Natural convection 
is further divided into cases of free convection and contained natural convec- 
tion. Natural convection implies that no causes of the motion exist other 
than effective body forces acting on portions of the fluid with density varia- 
tions due to thermal expansion. In that case, effects of pressure and appar- 
ent body forces can be expressed entirely in terms of buoyancy forces. The 
fluid convection associated with the present problem is generally a combina- 
tion of natural and forced convection. In the special (or limiting) case of 
steady rotation it reduces to the state of contained natural convection. 

Whenever a fluid flow is at low velocity relative to boundaries, with 
relatively small accompanying temperature and density variations, one usually 
considers the highly simplifying assumption of ’’constant density" of the fluid 
flow. A complication arises, however, when temperature differences, producing 
density variations that are acted upon by body forces, are a significant 
cause of the convection. Obviously, density variations must be accounted for 
in the body- force terms (the momentum production) , if not in other aspects of 
compressibility of the fluid flow. The usual approximate method of treating 
problems of nearly constant-density flow with "buoyancy" body forces (cf. 
ref. 1, p. 412; ref. 2, p. 320; ref. 3, p. 248) is commonly known as the 
Boussinesq approximation (e.g., ref. 4, p. 506; ref. 5 5 P* 684; ref. 6, p. 29) * 
and has been applied to natural convection driven by constant body forces. 

That approximation has recently been generalized to include spatially varying 
body forces in a uniformly rotating container (see, e.g., refs. 7-9)- 
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The Boussinesq approximation and its generalization to variable body- 
forces apply only to pure natural convection, and are not sufficiently general 
for present purposes (with time-varying 'rotation rates producing forced con- 
vection). Forced convection in rotating flows of homogeneous fluids has been 
studied extensively (e.g., refer to ref. 8), especially with reference to the 
"spinup" and "spindown" problems. The early fundamental treatments by von 
Karman and Cochran and by Bodewadt (see ref. 10, pp. 157-162) showed fundamen- 
tal differences between the spinup and spindown problems in three dimensions 
because of instabilities. A recent study by Briley and Walls (ref. 11 ) indi- 
cates that these effects also occur in a tank rotating about its center and 
lead to Taylor vortices in the spindown problem at some Reynolds numbers. The 
rotating forced- convection flows that have been studied did not include strat- 
ification effects, that is, body forces due to density variations. 

In the present study, a small-density- variation approximation is devel- 
oped in a manner analogous to the Boussinesq approximation, including both 
variable buoyancy body forces in a stratified fluid and the allowance for 
forced convection produced by time- varying rotation. Thus, the formulation is 
developed both to be consistent with the physics of the time-varying rotation 
(forced convection) and to reduce to a formulation for pure natural convection 
in the limiting case of steady rotation that is analogous to the Boussinesq 
approximation . 

The development procedure is to first write the compressible Navier- 
Stokes system of equations in vector form in an arbitrary noninertial frame of 
reference. Then the set of equations is reduced to an approximate form for 
nearly constant density flow, accounting properly for all body- force terms in 
the limit as the temperature and density variations vanish. (it will be seen 
that both the momentum and energy conservation equations require special con- 
sideration for the present problem.) For convenience in numerical analysis, 
the equations are then cast in terms of the vorticity and stream functions. 

The special case of a constant inertial axis of rotation and the further spec- 
ialization to two-dimensional flow are formulated. The two-dimensional equa- 
tions in rectangular coordinates are presented for use in computation of con- 
vection in a square tank. 

Although three-dimensional effects such as axial flow due to Ekman-laver 
suction (see, e.g., ref. 8) and development of Taylor vortices (e.g., ref. ll) 
cannot be represented in the two-dimensional formulation, the first numerical 
treatments of this problem in two dimensions are expected to yield significant 
information regarding the adequacy of vehicle rotation to produce sufficient 
convection. Because additional convection mechanisms are supplied by three- 
dimensional effects, the results from the two-dimensional formulation are 
believed to be conservative. 

COMPRESSIBLE NAVIER-STOKES EQUATIONS IN A TIME-VARYING 
ROTATING REFERENCE FRAME 


The Navier-Stokes set of equations for compressible flow in a noninertial 
(rotating and accelerating) reference frame is 
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DY 

Dt 

De 

Dt 


where D/Dt is the substantial-derivative operator 

— = — + V • v 

Dt at ~ - 


SB- + pv • V = 0 

(2.1) 

= -Vp + V • T + pg 

(2.2) 

= -pV • V - V • q + y$ 

( 2 . 3 a) 


and where 5 for later use, equation (2.3a) may also be written in an equivalent 
form, in terms of specific enthalpy (ref. 2, p. 322), as 


P 


Dh _ Dp 
Dt Dt ~ Y 


q + y$ 


(2.3b) 


The set of equations must be supplemented by appropriate equations of state: 

p = p(p,T) 

e = e(p,T) = / c v dT 
h = h(p,T) = / c p dT 

In the above equations, the Navier-Stokes expression for the viscous stress 
tensor t, defined by 

g = -pj + t (2.5) 

(where p is the thermodynamic pressure, I is the unit tensor, and g is 
the total stress tensor), is 



T = x(v • v)i + y[v Y + (vy) ] (2.6) 

u 

where 

y shear viscosity coefficient 

2 

A = tc - ^ p "second viscosity coefficient" 

k = ~^lA3 . )l , „ . : 2 & -bulk viscosity coefficient 

v • y 

(VV) "transpose" of the dyadic VV 

~ t 
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Also in equations (2.3), the viscous dissipation function 0 is defined by 

*= t : VV (2.7) 

and the heat-flux vector q is given by the Fourier heat-conduction law: 

q = -kVT (2.8) 


Of primary concern here 
noninertial system g, which 


is the apparent body force per unit mass in the 
is given by 


where 


g = -[a + ft x r* + 29 x v + n x (n x r«)] 


( 2 . 9 ) 


ft time-dependent angular velocity of the noninertial system rela- 

tive to the inertial system 

ft S dft/dt 

R# position vector from a fixed point (0* on fig. 2.1) on the axis 

of rotation of the noninertial system to a point of interest 
in the rotating system 


ft x (n x R*) 


centripetal acceleration due to rotation of the system 


2ft x v Coriolis acceleration 

ft x R* linear acceleration due to the angular acceleration ft 


a = a(t ) 


(time-dependent ) linear acceleration of the noninertial system, 
in excess of ft x R* 


As an effective body force per unit mass , a may be taken to include any body 
forces acting on the gas in the inertial reference frame. 

Equations (2.l) through (2.3) are derived directly from the corresponding 
equations in an inertial reference frame by substituting for the term 
p (DV/Dt ) , measured relative to the inertial system, the terms 



+ ft 


xR* + 2ft*y + ft>< (ftxR*) + a 


(cf. Becker, ref. 12, p. 251 )• The term p(DV/Dt) is the only term in the 
inertial system affected by the transformation to a noninertial reference frame 
because it is the only time derivative of a vector function (ref. 12 ). 

At this point it is convenient to define a position vector R as the per- 
pendicular vector from the axis of rotation to any point of interest determined 
by R* (see fig. 2.1), such that 



Q • R = 0 


( 2 . 10 ) 



Figure 2.1 Position vectors in a rotating system. 


The vector R always lies in a plane 
of rotation. Its origin 0 moves 
along the axis of - rotation as the point 
of interest inside the tank is varied. 
Note that 

n X R* = Q X R (2.11a) 

h x R* = ft x r (2.11b) 

where the axis of rotation is taken to 
be constant in inertial space. With 
equations (2.1l), R* may be replaced 
by R everywhere in equation (2.9) * or 


g = -[a + ? x R + 25 x Y + x (Q x R)] (2.12) 

REDUCTION OF EQUATIONS FOR SMALL-DENSITY-VARIATION APPROXIMATION 


Equations of State and Small-Perturbation Theory 

It is assumed that pressure and temperature variations , and therefore all 

fluid property variations, are small over the space of the tank and relatively 

slow in time. The thermal equation of state p = p(p,T) is then written as a 

Taylor series about the nearly constant average state defined by 

p = p(p ,T ) : 
o K *0 o 


= p + 

fe) 1 

AT + 

Vie.) "1 

M o 

\ ST / 

l_' / p - 

'0 

_l 3P / T _ 


Ap + 0[(AT) 2 ,(Ap) 2 ] 


= p Q {l - 3 AT + a Ap + 0[ ( AT) 2 , (Ap) 2 ] } 


where 


3 = 


=1 Ifi. 

P \3T 


P' 


is the thermal expansion coefficient at the state p Q , T , 


a = 


1 / 3£_\ 

Lp \ 3 p/ t J 


is the compressibility of the fluid at the state p^ , T q s and 


(2.13) 


(2.14a) 


(2 . l4b ) 
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(2.15a) 

(2.15b) 


AT = T - T 

o 

Ap = p - p Q 

As is customary* in convection problems involving buoyancy forces due to 
thermal expansion, when pressure gradients are expected to be very small, we 
neglect not only terms of 0[ ( AT) 2 , ( Ap) 2 ] in equation (2,13), but also neglect 
a Ap in comparison to g AT: 


p : p (l - B AT) (2.16) 

o 

From the resulting flow computations made using equation (2.16), one can 
always check that a Ap/g AT is small everywhere to justify £ posteriori the 
use of (2.16). 

By the same procedure as in equations (2.13) through (2.16), the caloric 
equations of state in (2.^) are simplified by assuming the pressure dependence 
of e and h to be negligible, so that the differentials de and dh can be 
written simply as 


de = c dT 


dh = c dT 
P 



(2.17a) 

(2.17b) 


The coefficients a, g, c , and c must be determined from thermodynamics 
and are not considered further V in the present article. 

The formulation involving the flow equations above with equation (2.16) 
substituted everywhere for p can be regarded as a small-perturbation theory 
with 

p = p (l + 6) (2.18a) 

or 

P - P_ 

6 = = ~g AT (2.18b) 

^o 

where 

1 <s ( « 1 (2.19) 

so that 6 can be neglected in comparison to unity. It will be important in 
formulating the theory not to neglect p 6 in those instances where after sub- 
stitution of (2.l8a) the term p Q effectively cancels out of the equations. 

Along with the general assumption of small temperature, pressure, and 
density variations, the transport properties y and k are taken to be 
constants at a given time. 
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If desired, the small-perturbation theory can be extended formally to 
higher approximations that supply higher-order corrections for significant AT 
and Ap in the Taylor expansions (about the state p Q , p Q , T Q ) of p, h, and 
the transport properties. The formal small-perturbation theory (with appro- 
priate emphasis on buoyancy effects to make them of lowest order) yields a 
first-order formulation equivalent to that derived below in a less formal way. 


Mass-Conservation and Momentum Equations 

In the limit as 6 0, one finds conservation of mass (eq. (2.1)), to be 

approximated by 

V • Y = 0 (2.20) 

This equation can henceforth be used in the flow calculation to represent con- 
servation of mass and can be combined with the other equations wherever 
appropriate . 

In deriving the approximation for the momentum- conservation equation as 
5-^0, one could proceed in a manner equivalent to that for the Boussinesq 
approximation only if the angular velocity of the tank were constant. In the 
present problem, however, we must allow for time-varying angular velocity. 

Thus, convection in a rotating tank, with fluid-density variations caused by 
temperature variations, falls essentially into two categories: 

1. If the angular velocity of the tank is constant 3 then the conditions 
are met for the limiting state referred to as "contained natural convection." 
For that case, in the limit as 6 -* 0 , the fluid approaches a state in equilib- 
rium with the walls (no relative motion), which is equivalent to a rigid body 
rotation. If, then, there are small temperature and density variations in the 
tank, the resulting flow is entirely a small perturbation on the "rigid-body" 
rotation. In that case, no causes of the motion exist (in the reference frame 
of the tank) other than the effective body forces acting on portions of the 
fluid with density differences due to thermal expansion. Then the effects of 
pressure and apparent body forces can be expressed entirely in terms of buoy- 
ancy forces, with essentially constant density of the fluid assumed in all 
other respects. These are the assumptions of the Boussinesq approximation used 
frequently in special cases of natural convection. 

2. If only the axis of the angular velocity is constant (in inertial 
space), but the magnitude of rotation changes significantly with time (e.g., a 
steady rotation is stopped or reversed) , then there are factors producing the 
convection other than simply the apparent body forces acting on density dif- 
ferences due to thermal expansion. Sudden changes In the boundary rotation 
from a previously steady rotation produce boundary layers, with significant 
velocities at the boundary relative to the internal core. The rotational 
inertia of the core tends to keep it rotating without substantial change, 
except for factors such as buoyancy body forces and slight axial flow due to 
Ekman-layer suction (e.g., see ref. 8), until affected by diffusion of the vis- 
cous forces from the boundary. In this case, the assumptions of natural 
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convection are violated, and a more general formulation must be used. It 
should, however, for present purposes, include the valid approximate formula- 
tion of natural convection in the special case where the rotational velocity is 
held constant for a period of time. 

In the limit as 6 -> 0 (p ~ p Q ) , y and k approach constant values at a 
given time, and from (2.6) and (2.20), 

V • x = y[V • vy + V • (VV) ] 

= y[V 2 V + V(y • V)] * yV 2 V (2.21) 


Assuming that only the apparent-body- force term in (2.2) may be affected by 
small density variations, we therefore retain 6 only in that term and write 
(2.2) as 

DV i , 

= " — YP + ^ Y + (1 + «)g (2.22) 

Po 

By manipulating equation (2.22) in various ways, one can determine in what 
sense 6 must he retained in (2.22) in the limit as 6 0. For this, various 

terms in the quantity g, given by equation (2.12), can be expressed as gradi- 
ents of scalar quantities. Those portions of (l + fi)g that can be so 
expressed can then be regarded simply as modifications to the pressure gradient 
term, (-l/p Q )Vp (which is ultimately to be eliminated from the problem). The 
various forms equivalent to (2.22) are derived in this and following sections. 

First consider the centripetal acceleration, Q x (Q x R). For arbitrary 
vectors A and B the following vector identity holds (ref. 13* p. 270): 


A x (B x C) = (A • C)B - (A • B)C 


Therefore 


Q x (ft x R) = (ft • R)Q - 


(2.23) 

(2.2U) 


One can write 

£2 = e3fi(t) 


(2.25) 


where §3 is a unit vector in the direction of the axis of rotation, so that 
with ( 2 . 10 ), equation {2.2k) becomes 

n x (n x r) = _n 2 R = -y[(i/ 2 )^ 2 r 2 ] ( 2 . 26 ) 

Further, consider the fact that since a is , at most, a function of time, 
a = a(t ) , then 

Y x a(t) = 0 (2.27a) 
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so that a scalar potential <f> = R* % a(t) can be defined such that 

a = V<J> (2.2Tb) 

satisfies the vector identity 

V x v<f> = 0 (2.27c) 

Then equation (2.22), with (2.12), (2.26), and (2.27b), becomes 

DV i 0 , . » 

r-r = - — VPi + vV 2 V + g6 - (ft x R + 2ft x v) 

Dt P 0 ~ - ~ 

where 

Pi = p + P Q [<f> - (l/2)ft 2 R 2 ] 

(with and ft both depending on time), which may be called the "reduced 
pressure" (cf. ref. 8, p. 6, for constant ft). 

On the Use of Stream Functions to Reduce the Coriolis Term 

Under some conditions, the following procedure introduces further 
simplification: 

For any arbitrary scalar functions and ^2 th e vector identity holds 

(ref. 13 , p. 278) : 

V • (Vi/jj x v^ 2 ) = 0 ( 2 . 30 ) 

To satisfy (2-30), a vector V* (to be defined for special cases) can be writ- 
ten in terms of ^1 and. ^2 as 

V* = x V\p 2 (2.31) 

(The vector Y* is normal to both and V^ 2 ? and is therefore along the 

intersection of the two surfaces = Ci(t) and \p 2 = C 2 (t). In cases where 

one defines Y* to be V, equation (2.20) is satisfied by (2. 30); the surface 
intersection is then a streamline and ipi and \p 2 are "stream functions." 

With equation (2.3l) and use of the identity (2.23), the Coriolis accel- 
eration can be written as 

2ft x Y = 2ft x (V\p 1 x V^ 2 ) + 2ft x (V - V*) 

= 2 [ ( ft • 7jp 2 )Vifj 1 - (ft • V^ 1 )Vi^ 2 ] + 2ft x (Y - Y # ) (2.32) 

Use of this expression is convenient only in some special cases. 


(2.28) 

(2.29) 
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Case 1 . No special restrictions — With no special restrictions, it is just 
as well* to use ( 2 . 28 ) without (2.32) because use of ipi and ^2 (2.32) 

appears to offer no simplification. 

One considers then in (2.28) (.in the limit as 5 0) omitting the parts 

of g in the term gS that are the same as the remaining body forces: 

-h x R - 2ft x V. However, if it should happen (as is true in some cases) that 
the Coriolis force term (-20 x V) is exactly balanced by part of the pressure 
gradient and is effectively eliminated from the problem, then that part of g 
should not be omitted from the term g5 for small 6. Therefore, if the 
effects of 20 x v are unknown, it is better to leave 20 x V in g in the term 
of order 6, as well as in the body force term without the factor S. 


Case 2 . No flow changes in z direction — Let 63 te a unit vector in the 
z direction, also the direction of 0(t) (eq. (2.25))- Call the plane z = 0 
the xy plane (although other coordinates may be used in the plane). Let V* 
be the velocity projection in the xy plane, so that 

V = Y*(x,y,t) + e 3 w(x,y,t) (2.33a) 

Note that V * V = V ♦ Y* = 0- Then let 


4»2 = z > ih = iKx,y,t) 
so that (for use in (2.32)): 

V^2 = §3, £} • Vijj2 = fi(t) 

In rectangular coordinates (x,y,z): 


Y<h 



§2 


3y 


V* = Vipi X V\p2 = ej 


9y 



(2.33b) 


(2.33c) 


(2.33d) 


But, by definition, V* = eju + e 2 v so that if the flow is two-dimensional, 4> 

is the conventional stream function. In cylindrical coordinates (r,0,z) with 

unit vectors e , e_, e?: 

~r ~ 0 ~ 0 


Vi|/ 1 


e ? 

~r 9r ~0 99 




Y * = vii > 1 x v\jj 2 


l dih 
~r r 90 


e 

~6 9r 


> (2.33e) 


= e v + e„v„ 
~r r ~0 0 


J 
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Also 


2 • V\J»! = ° ^ 

Y - V* = e 3 w and 2 * (V - Y*) = 0 J 

> (2. 33f ) 

With equations (2.33), equation (2.32) "becomes 


22 x Y = V[22(t)ij>] 

(2.3U) 

and ( 2 . 28 ) becomes 


DV 1 , . . , 

~ = - — VP 2 + vV 2 V + g6 - 2 x R) 
Dt P o~ ~ ~ 

(2.35) 

where 


P 2 = P + P [<t> - (l/2)2 2 R 2 + 22t|/] 

(2.36) 


with <f> and 0, time dependent; see reference l4, p. 177 * for <J> and 
independent of time . 

For the case where the rotational speed is constant (ft = 0) , equation 
(2.35) is analogous to the Boussinesq approximation, with pure natural convec- 
tion caused "by M buoyanc y" body forces. 


The case where the z component of velocity w(x,y,t) is zero {two- 
dimensional flow) is included in case 2 , equations ( 2 . 33 ) through ( 2 . 36 ). 


Case 3. Rotational symmetry 1 in cylindrical coordinates — Consider cylin- 
drical coordinates (r, 0 ,z), with the unit vector e„ = e^ = ft/ft defining the 


z direction and the coordinates r and 


z 



Figure 2.2- Cylindrical coordinates 
for rotational symmetry. 


0 being in the xy plane (fig. 2 . 2 ). 
Let e r and §q be unit vectors in the 
respective directions of increasing r 
and 0. Assume that the flow is 
invariant in 0. (This would be true 
in the rotating tank only if the axis 
of rotation is at the center of the 
tank and only for those ranges of the 
parameters for which such a flow is 
stable . ) 

Let any plane through the axis 
(constant 0 ) be called the rz plane. 
In this case, let V* be the projec- 
tion in the rz plane of the velocity 

V at a point in that plane, so that 

V = V*(r,z,t) + e 0 v Q (r,z,t) (2.37) 


*A distinction is made between rotational and axial symmetry (v 0 = 0 ) fol- 
lowing Synge (ref. 15) and reference l 6 . 


21 



Note that with rotational symmetry 


V • y = V 


F = ;e (rv ) 
r dpF r 



( 2 . 38 ) 


where v r and v z are the respective velocity components in the r and z 
directions, and do not depend on 8 . 


Then let 


i ^2 = 9 and = ^(r,z,t) 


so that, for use in ( 2 . 32 ) 


1 S *2 1 


Viio = e rt — = — e 


0 r 86 r ~0 


0 * YF 2 = 0 


V*i 

V* 


dip dip 

f ^ i ■ -U a 

~r 8 r ~z 8 z 


= x ^2 = e r ^- + § z (f If 


= e v + e v 
~r r ~z z 


ft * = -ftrv 

Q x (V-V*) 


-e 

~r 0 


A 




(2.39a) 


(2.39b) 


Therefore, from (2.32), with (2.39b), 


2f2 x v = 2ft (v e Q - v Q e ) (2.U0) 

" - r~t) y~r 

It Is seen from this special case of rotational symmetry that the Coriolis 
force 2ft * V in equation (2.28) is in general not negligible in three- 
dimensional flow, either when 6 -* 0 with forced convection (ft ^ 0 at some time 

in the problem) or for steady rotation (ft = 0) with 6^0. Thus, if the flow 

is restricted to be two-dimensional, then (see case 2 above) the Coriolis force 

can be included in Vp and eliminated in the limit as 6 -* 0 ; but if the flow 

is three-dimensional, the Coriolis force is not balanced by pressure gradients, 
and so cannot in general be eliminated when 6 0 . 

Case 4 . Rotational symmetry 2 in spherical coordinates — Following a proce- 
dure similar to that in case 3 above, consider spherical coordinates (r,(j), 0 ), 
with the unit vector §3 = ft/ft defining the axis of symmetry from which the 


2 See footnote 1 . 
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Figure 2.3.— Spherical coordinates for 
rotational symmetry. 

Note that 


polar angle <p is measured (fig. 2.3). 
The assumption of rotational symmetry 
implies the flow is invariant in the 
azimuthal coordinate 0, although there 
may “be a velocity component 
v 0 =. v 0 (r,<j>,t) . 

Let e r , e ( p 9 and §q he the unit 
vectors in the respective directions of 
increasing r, $ , 0. Let any plane 
through the axis (containing e r and 
e^, as well as §3) he called the r<p 
plane. Let V* "be the Projection of 
Y(r,4>,t) at the point (r,<f),0) onto the 
r 4> plane, so that 

Y = V*(r,(j),t) + e 0 v 0 (r 

(2.41) 


V = V • V* = g7[r 2 v r (r,ij),t)] 


r — 7 tx [(sin <j>)v (r,<|>,t)] = 0 

r sin <|> 3<(> 4> 


(2.42) 


where v r and v^ are the respective velocity components in the r and <f> 
directions, and do not depend on 0. Then, for use in (2.32), let 


so that 


From (2.U3h), 


and so 


Noting that: 


ip 2 = Q and = ifj(r,<J>,t) 

- I 1 -)e 

\ r sin 4>/ ~ 0 


V\p 2 

U • V\p2 = 0 

Vih = e ?• + e - 
1 ~r 3r ~<p r dtp 

V . V*, X V * 2 = If) * e t ( 


m 


r sin 3r; 


v = 


= e v ( r ,<J> ,t ) + e v (r 9 <J>,t) 
~r r ~<f> <f> 


-1 3± 




1 It A 

r r2 sin <p 3(p an % r sin <j> dr 


Vipi = e (-v.r sin <j>) + e,(v r sin ({>) 
~ ~r cp ~<p r 


§3 % = cos ^ 

es ' e = cos[(ir/2) + 4> ] = -sin <J> 
~ (P 


(2.43a) 


\ (2.43b) 


(2.43c) 
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we have 


(2.43d) 


Also vith 


R 


V^i = -R(r sin ij>)(v^ cos + sin <j> ) 
§3 x §0 = -(§ r s i n $ + cos (jj) 


we have 

fl x (V - V*) = -Rv„(e sin <b + e, cos 4>) (2.1+3e) 

~ 0 ~r ~<p 

With equations (2.43b, d, and e), equation (2.32) gives then 

2Q x Y = 2ft[-v 0 (sin 4>)§ r - v 0 (cos <{>)e^ + (v^ sin cj> + v^ cos 4> ) e ^ ] 

9 ~ v ; (2.44) 

For flow vith rotational symmetry, equation (2.44) for spherical coordinates is 
equivalent to equation (2.40) for cylindrical coordinates. Comments on case 3 
also apply here. 

Further comments on the proper forms of the equations to use for small 6 
can be made when the respective forms of the vorticity equation are derived for 
the various cases in the next section. 


Vorticity Equation and Vorticity-Stream Function Relations 

Computation of a flow field with use of vorticity as a primary variable is 
advantageous because (see Lighthill, ref. 17, pp. 57-60; and Greenspan, ref. 8, 

pp. 20-21): 

1. The unknown pressure and all conservative body forces are eliminated 
from the problem; 

2. Large, sudden changes in velocity or angular velocity of the surface 
produce large sudden changes in fluid velocity and large impulsive pressures , 
whereas the vorticity distribution varies smoothly. (Vorticity changes are not 
propagated at the speed of sound, as are velocity and pressure changes.) 

The fluid vorticity u) is defined by 

u) = V x v (2.45) 

To obtain the vorticity equation corresponding to each of the momentum- 
equation forms (2.28) and (2.35) use first the vector identity 

v • yy = v[(i/ 2 )v 2 ] - v x (v x y) (2.46) 

in the left side of each of those equations, and then take the curl of each 
equation. With use of several identities and with g given by (2.12), one 
obtains from ( 2 . 28 ): 


ii 
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9w 

at 


r - V x (y x (j) - vV 2 y = V x (-fi x R - 2S2 x V + g6) 


(2.47a) 


and from (2.35) for flow independent of z: 


3oa 


77* - V X (y X co) - vV 2 y = V X (-Q X R + g6 ) 


3t 

These forms of the vorticity equation are reduced further as follows: 
vector identities 


(2.47b) 
With the 


and 


7 x (A x B) = B • VA - A • VB + A(v • B) - B(V • A) 
V*aj = V*(VxV)=0 


(2.48) 

(2.49) 


and with use of the approximate mass-conservation equation (2.20), the left 
sides of equations (2.47) become 


„ Du , 

tt- + V • Vw - w • VV - vV 2 io = ^rf" - a) • VV - v 7 z oj 
dt ~ ~ JJt 


(2.50a) 


(The significance of u> • VV is discussed by Batchelor, ref. 14, pp. 267-268. ) 
It may be noted for later convenience that with use of the following vector 
identity for the dyadic AB, 

V • (AB) = A • VB + (V • A)B 

along with (2.48), one can also write the left sides of equations (2.47) in the 
"divergence form" 


TT + V • ( Vw) - V • (o)V) - VV 2 0) 
dt ~ ~ ~ ~ 

where, for example, in Cartesian coordinates 

~ * = 9x + 8y + 3z" 

? • <»y> - £ <« x y) + ^ '“yY* + si ( “zY> 

The right side of (2.47a) is reduced further by using the identity 

V x (g6) = 6 (V x g) + (76) x g 
with equations (2.12), (2.2 6), and (2.27a) so that 


(2.50b) 


(2.51) 
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7 x [-fl x R - 2fl x V + giS] = (l + 6)V x (-fl x R - 2g x V) + (VS) x g 


(2.52) 


If then the quantities V x (2 x R) and V x (ft x y) are expanded according to 
(2.48), for the case under consideration (constant inertial axis of rotation) 
•with 


2 = e3^(t) , V2 = 0 , V • 2 = 0 
2 • VR = 0 , V • R = 2 


(2.53) 


equation (2.47a) with (2.50a) and (2.52) becomes 


where 


D(jO 

- v * YY ~ vv2 “ = (Y 6 ) x g + 2(1 + s )(2 • vy - 2 ) 


ft • VV = 2(3V/3z) and 2 = <532 


(2.54) 


We have now arrived at a point where, for consistency with the previous develop- 
ment for small 6 , we can neglect 6 in comparison to unity in the body- force 
term. Equation (2.54) is then approximated finally by 


(VS) x g + 2(2 • VV - 2) 

(VS) x g + 2(ft3V/3z - e 3 2) 

In the case where the flow is -independent of s 3 equation (2.55) 
to equation (2.47b) reduced to 

57 " ~ y * YY - = (Y 6 ) x g - 53 ( 22 ) 


Du 

Dt 


- to 


VV - vV 2 io = 


(2.55) 
is equivalent 

(2.56) 



In the further specialization to tuo-dimensionaZ fZow (w = 0), the vorticity is 

to = §3(0 (2.57) 

so that (2.56) becomes the scalar equation 


— _ vV 2 to = [VS x g ] _ 22 
Dt ~ 17 


~ z 


(2.58a) 


or, with the left side replaced by (2.50b), in two-dimensional Cartesian 
coordinates (x,y). 
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3 03 9 r \ 9 / \ n 9 96 96 

_ + _ (uu ) + — (VU)) - W 2 0> = g 2 ^ - gl — - 


(2.58b) 


or in two-dimensional polar coordinates (r,^), 


3(0 

at 



(rv oi) 
r 


+ 9? (v e w) 


vV 2 o) = g 


9_6_ 

0 9r 



(2.58c) 


where . gi and g 2 are the x and y components of g, and g r and g^ are the 
r and 0 components of g in (2.12). 


It is noted from the vorticity equation (2.55) that: 


1. In the limit as 6+0, the nonsteady rotation term in the momentum 
equation ( 2 . 28 ) is exhibited only in the term -2ft in the vorticity equation, 
whereas the Coriolis term results only in the term 2ft • VV in the vorticity 
equation, and both of these are in general not negligible as 6 + 0 in three- 
dimensional flow where 9V/9z ^ 0. 

2. In the limiting case of rigid-body rotation (relative velocity V + 0 
and 6+0), equation (2.55) reduces correctly to 



3. In the special case of steady rotation (ft = 0 for all time), the only 
driving term in (2.55) or (2.56) is the buoyancy-force term (V6) x g. 

4. In general, the driving factors in (2.55) are both the buoyancy term 
(V6) x g and -2ft, these being the only terms that are independent of, or con- 
tain terms independent of, the relative velocity y. Thus, if both 

(V6) x g = 0 and ft = 0 , we have V + 0 so that a rigid-body rotation (condition 
of no convection) is approached. The Coriolis term 2ft • VV in the three- 
dimensional equation plays an essentially passive role; that is, it apparently 
acts simply to absorb (oppose) part of the effect of the driving terms (see 
Ostrach, ref. 7), since it contains V and. vanishes as the flow approaches a 
rigid-body rotation when both (V6) x g and ft are zero. 

In connection with the use of the vorticity equation (2.55) or (2.56) it 
is convenient to use the stream functions for the special cases noted in the 
previous section, and to relate the stream function to the appropriate 
vorticity component. 

In case 2 above, with no flow changes in the z direction, y, defined by 
( 2 . 45)9 becomes 

U) = V X [e 3 w(x,y ,t ) ] + V x {v x [e 3 ij>(x,y ,t) ]} 

From the vector identity 

V X (v X a) = v(v • A) - V 2 A 

one has 

V X [V X (e 3 i|>)l = -e 3 V 2 iJj (2.6l) 


(2.59) 

( 2 . 60 ) 
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so that in rectangular Cartesian coordinates 


3w 3w „2 1 

y = SI e 2 e 3 V z * 


In the case of cylindrical coordinates (r,0,z) with flow independent of z, 


/I 9w\ 9w _2 r 

9 = W 5?) - h 37 - «3V 2 * 


where 


V2. = _i r _i L m + ifi 
V v r2 r 9r r 8r 302 


For two-dimensional flow, with (2.57) for w = 0, (2.62) and (2.63) reduce to 

a) = -V 2 * (2.6^ 

In case 3 above, with rotational symmetry in cylindrical coordinates, 

Y = e v ( r ,z,t) + e v (r,z,t) + e 3 v (r,z,t) (2.6^ 


where 


03 = V x V = e 0) + e -w. + e 3 oi 

~r r ~0 0 ~ J 2 


o) = - “v— 
r 9z 


9v 9v 
z^ r 

9r 9z 


19/ \ 

03 = — — ( rv n ) 

z r 9r 0 


19 * , 13* 

v = - — and v = — 
r r 9z z r 9r 


from (2.39b);, equation (2.66c) becomes 


JL I 11 ] , Ill 

r 9r r 9r / 9z 2 


3r2 r 9r 9 z 2 


Note that in this case, w is not directly related to V 2 *, since 

0 
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(2.69) 


V 2. = I _i / r 3 l\ + ifl 1 

^ r 3r 1 dr) 3z2 I 

= l£i + i lt+ ifl f 

3r2 r 3r 3z2 J 

Equation (2.68) would be useful in conjunction with solving (2.55) for this 
special case. 

In case 4 above, with rotational symmetry in spherical coordinates 


where 


Y = S r v r (r,<|>,t) + S ( r * <f> »t ) + e 0 v 0 (r 


» = ®r u r + V+ + §6 w e 


(2. TO) 
(2.71a) 


With 


w = ~ 1 -V ~ . ' TT~ [(sin <f>)v ] 
r r sm q> dtp 0 


% = "T sf (rv e } 


-1 

r 


'Hr : 

9(j> 9r 


i M d = -i 
v r r2 sin <J> 9<j> 311 V <f> r sin <f> 9r 


v = 


(2.71b) 

(2.71c) 

(2.714) 

( 2 . 72 ) 


from (2.43c), equation (2.71d) becomes 


-r 3 (sin *)(o 0 = r 2 + (sin 4> ) ^ 


1 11 

sin (f> 9<j) 


(2.73) 


Note again that w is not directly related to V 2 i)j, which is given by 

0 


r 2 v 2^ = _L ( r 2 ll\ + _1 1 f( sin + ) 

v 9r \ 9r j sin <f> 9cJ> r 91 9<f> 


(2.74) 


Equation (2.73) would be useful in conjunction with solving (2.55) for this 
special case. 


Energy Equation 

The energy-conservation equation is considered in two forms, (2.3a) and 
(2.3b). With the caloric equations simplified to the forms of (2.17) with 
V • Y assumed to be zero (eq. (2.20)) in expressing the viscous stress tensor 
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t in (2.6) for use in (2.7) so that 

4> = [TV +.(VV) t ] : VV (2.75) 

and with k in (2.8) assumed to be nearly constant as noted above, the energy 
equations (2.3) take the forms 

pc ir + PV • V = kV 2 T + (2.76a) 

v Dt ~ 

and 

pc 22. - = kV 2 T + y$ (2.76b) 

p Ut DZ 


For the conditions under consideration in this problem, we have retained the 
term pV * V in (2.76a) for reasons given below. 

If there is any apparent conflict in the results from either of the two 
forms of the energy equation, one should use the form that is most consistent 
with (2.16), rather than with the more restrictive equation (2.20), which 
results from letting 6 -*• 0. Obviously, for a liquid in which Cp ~ c v and 
p z constant, there is no difference in the two equations (with Dp/Dt very 
small). But in the present problem, there is a significant difference between 
Cp and c v , so we must consider the relative consistency of either neglecting 
pV ' V in (2.76a) or neglecting Dp/Dt in (2.76b). 

Consider first equation (2.76a). One could use (2.76a) with (2.20) sub- 
stituted if, in fact, 

|pV Yl « |kV 2 T | 

To determine whether pV * V is negligible, write (2.l) as 


pV • V = 

and use the approximation from equation 


_ J2. 2ji 

p Dt 

(2.16) that 


1 Dp -B DT DT 

p Dt 1 - g AT Dt * B Dt 


so that 


DT 

pV • V = P B — 


(2.77) 


Near the critical point of a fluid, for example oxygen with T * 150° K and 
6 - 0.01/°K, one finds, for example, from the van der Waals equation of state, 
that 


P 

P 


> 


(RT) 


perfect gas 


= 0(c v T) 


( 2 . 78 ) 


k 


j & 
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Therefore, from equations ( 2 . 77 ) and ( 2 . 78 ), 


or 


Ip? 


(BT)pRf 


pc 


DT 
v Dt 


Ip? 


pc 


DT 

Dt 


( 2 . 79 ) 


The result ( 2 . 79 ) ' indicates that pV • V would not be negligible in determin- 
ing DT/Dt from (2.76a). 

Consider then the use of equation (2.76b), with use of the approximate 
form of the equation of state (2.16), where (3 (eq. ( 2 .l 4 a)) is to be deter- 
mined from a functional relationship such as 


p = p(p,T) 


Because of ( 2 . 80 ), we can write 


5R - 

Dt 



With the approximation represented by 


2e. + t*s\ 

Dt \ 9 T / Dt 
equation (2.l6) 


( 2 . 80 ) 


(2.81) 


2e. = 

Dt 


fl DT 
P o S Dt 


(2.82) 


and from the definition of 


g in ( 2 .l 4 a) , 


equation (2.8l) hecomes 



(2.83a) 


But for a state function such as ( 2 . 80 ), a fundamental identity from calculus 
states that the bracketed quantity in equation (2.83a) is identically zero. 
Therefore 


. 

Dt 


0 


( 2 . 83 b) 


This result shows in a rational way that use of the approximate state relation 
(2.16) is most compatible with the approximate energy equation (2.76b) in the 
form (as 6 -* 0 in (2.18a)): 


DT 

Dt 


(2.810 



(See also ref. 10, pp. 126, 127.) With use of (2.20), the left side of (2.84) 
may also he written as 


M = 3T () 

Dt 3t ~ ' 

where in the special case of two-dimensional flow, 


(YD = -jf (uT) + ^ (vT) 

= r 3? (rT r T) + r 36 U e T) 


(2.85a) 


(2.85b) 

(2.85c) 


PROBLEM DEFINITION AND APPROXIMATE FLOW EQUATIONS 
FOR TWO-DIMENSIONAL SQUARE TANK 


In this section, the theory developed above is summarized and specialized 
for computation of the convection caused by time-dependent rotation and density 
variations due to temperature variations in a two-dimensional square tank. The 
configuration is as shown in figure 2.4. The x,y coordinate system is fixed 
relative to the tank, with the origin at the corner as shown. The rotation is 
about a point a distance R c from the tank center, and the radius vector from 
that point is denoted by R (see fig. 2.1) with components in the x and y 
directions equal to 


Rl = x - \ 1 

R2 = y - \ £ + R c 



J 


( 2 . 86 ) 


The time-dependent angular velocity Q 
is positive when counterclockwise. In 
the x,y coordinates the effective body 
force per unit mass, given by equation 
(2.12), becomes 


g = 

§181 + 

§2g2 

(2.87a) 

gl = R^ 2 

+ 2Qv + 

R2^ - 

• a l 



(2.87b) 

g 2 = R2^ 2 

- 2Qu - 

R^Q - 

- Si 2 




(2.87c) 



Figure 2.4.— Square tank configuration. 
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and where &i and a .2 are prescribed x and y components of a (defined follow- 
ing equation (2.9)) and u and v are the x and y components of V (eqs. 
(2.33)). 

With 6 given by equation (2.18b), the vorticity equation (2.58b) becomes 


+ — — (uoi) + “■ ( vco ) = \> 
9t 9x v ' 9y v ' 


9 2 oj 9 2 w\ 0 9T . 9T 

9x2 + 3y2j g26 9x Sl6 9y 


- 2£2 


( 2 . 88 ) 


The energy equation (2.8^) with (2.85) "becomes 

9T + _9 (uT) + _i (vl) . _Jt./l£T 

at 3x l ' 8y ' P o C p \ 3x2 3y2/ Cp 

where the dissipation function, from (2.75) becomes 


(2.89a) 


$ = 2 


9u 
L\ 3x / 


r * (g 


W J 


+ i|i+ 

9x 9y 


(2.89b) 


It is convenient to use the vorticity-stream function relation from equation 

(2.6b): 


+ A = _ w 

9x2 + 9y2 


where from (2.33d), 


u = ^ 

ay * 


v = _ 14 

9x 


( 2 . 90 ) 


(2.91) 


Sufficient initial and boundary conditions are needed to specify the prob- 
lem completely. Initially (t = 0), a temperature distribution is specified and 
the velocities u and v are everywhere zero (rigid-body rotation). At all 
times , 


on x = 0 


y = 0 



v = 0 


or 


14 = 14 = 0 

9y 9x 


(2.92) 


The condition in (2.92) on the normal component of velocity (tangential deriva- 
tive of ip) can be replaced by ip = 0 on the boundary for application .to equa- 
tion (2.90)- The condition on the tangential component of velocity (normal 
derivative of ^) can then be incorporated into a condition on w at each 
boundary for equation (2.88). The temperature or its normal derivative may be 
specified on the boundary at all times for equation (2.89a). 

The function £2(t) must be prescribed as a condition of the problem. As 
one example, for sudden reversal of rotation (with £2 = -£2 0 for t < 0), one may 
use an approximation of 
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«(t) = -fi + 2Q S(t) 
o o 

where Q is a constant and S(t) is the unit step function, 


(2.93a) 


Then from (2.93a) 


S(t) = 0 , t < 0 

= 1 , t > 0 

ft(t) = 2Q S(t) 


(2.93b) 


(2.93c) 


where S(t) is the Dirac delta function. The simplest numerical approximations 
to S(t ) and S(t) may be used; for example (with very small). 


S(t) = 0 


for t < 0 


= t/ti for 0 < t < t] 


= 1 
S(t) = 0 


for t > ti 
for t < 0 


= 1/tj for 0 < t < tj 
= 0 for t > 1 1 


(2.9^a) 


(2.9^b) 


For convenience in further treatment, equations (2.86) through (2-9*0 may 
be put into dimensionless form using the following dimensionless variables and 
parameters (where L is an arbitrary length): 


£ = - 
* L 5 


n = 


_ £ 


T = 


vt 


* tj _ 

a)* = , H - 


L ’ 1 " L2 

T - T 0 

T r - T ’ 
n o 


V = 


_ i 


v 


U = — , 
v 


V = 


vL 

v 


D*(t) = 


2fl(t) 

vL-2 


2G 
= ~T~2 


(2.95a) 


Gl = Gr(gi/G^L ) , G 2 - Gr(g 2 /^L) 
a* = a^/^L , a* = a 2 /fi^L 


4>* = (L 4 /v 2 )$ 


J 


with 


3 *+ 



Reynolds number 


Prandtl number 
Grashof number 
Eckert number 


Re 

Pr 

Gr 

Ec 


ft L 2 
o 


p VC' 

-P- R 


(n2 L )L 3 6(T_ 

Q n_ 

V2 

(ft L) 2 
o 

c (T p - T ) 

P X\ O 


T o } 


\ (2.95b) 


J 


Equations (2.86) through (2.9*0 become, in dimensionless form (dropping the 
asterisks from co*, ft*, $*, a*, and a* from here on): 


3u> 3 / IT \ 3 /, r % 3 2 (d 

ji (U “ ) 8n (V “ ) = 352 


3r^ 




3H 3 , TTH x 3 ,. m x 1 /s 2 H . 3 2 H^ 

37 + 3? (UH) + 3^ (VH) ° P? ( 3i2 + ^2, 


_ Ec . 
+ Ri 2 $ 


9 2 ,j, a 2 ^ 

352 9 n 2 “ 


(2.96) 

(2.97) 

( 2 . 98 ) 


where 


U = — 
3h ’ 


V = - tt 


3? 


Gl = 


Gr 


Re 2 


5 - 


IrHft 2 + w + 


n - 


2 L L 


1 dft 

2 dx 


G2 " Re 2 


1 l ^ 
n ' 2 L + 


$ = 2 


with conditions 


3U 

3C 


3V 

3n 


/ 3V 3U\ 2 

^3£ 3n/ 


5 = 0 and £ = £/L 0 JW = 3Y _ ^ = q 

n = 0 and n = &/L J ^ ^ 


(2.99) 


- Gr aj (2.100a) 


- Gr a, (2.100b) 


(2.101) 


( 2 . 102 ) 


The value of H or its normal derivative may be specified on the boundary. 
The dimensionless function ft(x) may be specified, for example, by equations 
(2.93) and (2.9*+) with t and t^ replaced by x and xj. 
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CONCLUDING REMARKS 


A set of general approximate equations has been developed that represents 
the Navier-Stokes description of convection of a thermally stratified fluid in 
a container with arbitrary time-dependent rotation. The equations are valid 
for combined forced and natural convection with significant, but sufficiently 
small, density and temperature gradients. All relevant terms representing 
effects of rotation and changes in rotation are included. 

The special case of convection in a two-dimensional square container was 
formulated for subsequent computations . 
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3. COMPUTATIONAL METHOD FOE CALCULATING 
CONVECTION IN A ROTATING TANK 

Harvard Lomax and F. R. Bailey 



SUMMARY 


A time-dependent finite-difference method is presented for calculating the 
convection of vorticity and energy in a stratified, supercritical cryogenic 
fluid contained within a two-dimensional square tank rotating, in a gravitation- 
less field. The finite-difference approximations to the convective part of the 
governing small- density -variation form of the Navier-Stokes equations are "based 
on an explicit predictor-corrector scheme of second-order accuracy in time and 
space. A discussion of the stability and accuracy of this method is included. 
The solution for the stream function that appears in the governing equation is 
determined by using a direct solution of Poisson* s equation based on double 
cyclic reduction. 


INTRODUCTION 


In this chapter, a numerical finite-differencing scheme is described that 
can be used to compute the convection of vorticity and energy within a two- 
dimensional square tank rotating in a gravitationless field. 

The purpose of the analysis is to determine whether local temperature 
stratifications in an insulated oxygen tank traveling in space at high pressure 
could be broken up, or stirred, by very low rates of rotation and by rotation 
reversals. Effects of both viscosity and heat conduction are included in the 
study of the time-dependent mixing. The mathematical model that describes the 
two-dimensional behavior of such an environment is given by the Navier-Stokes 
equations for small-density variations, one form of which is given in the next 
section. Many numerical calculations of these equations have been carried out 
but not with the physical conditions described above. These physical condi- 
tions guided the choice of the numerical procedure. In particular, a method 
was chosen that has very little numerical dissipation but is not highly accu- 
rate with regard to dispersion, since the amount but not the details of the 
mixing was considered to be of paramount importance. 


PROCEDURE 


The geometry of the problem is shown in figure 3.1. A two-dimensional 
fluid within a square is being rotated at a rate ft about some point exterior 

to the boundaries of the square. 

Since we are presenting only the 
basic principles of the numerical 
process , the details of deriving the 
particular form of the Navier-Stokes 
equation given next are omitted. The 
governing equations derived in chap- 
ter 2 are 
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(a) Geometry 


e=o € =10 

(b) Finite difference grid 


Figure 3.1.— Geometry and differencing grid for 
two-dimensional square tank. 
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The independent variables t and (5,n) are dimensionless forms representing 
time and Cartesian (x,y) space, respectively; u>, H, and ¥ are the dependent 
variables representing dimensionless forms of the vorticity, temperature, and 
stream function, respectively. The variables Gi and G 2 represent the apparent 
body forces in the rotating system as defined in chapter 2 and Pr is the 
Prandtl number. Equations (3.1a), (3.1b), and (3.1c) are referred to as the 
vorticity, energy, and Poisson equations, respectively. Note that the vortic- 
ity and energy equations (later referred to as the transport equations) are 
written in conservative form. 


The finite-difference computational domain is shown in figure 3.1(b). The 
crosses represent the boundary points corresponding to the walls of the con- 
tainer, and the open points refer to the interior points at which the dependent 
variables can change values. Any grid point at which a dependent variable can 
change in time is referred to here as a moving point for that dependent vari- 
able. This simplifies the description of the matrix formulations used in the 
analysis of the numerical methods. The value of any dependent variable, say 
a), at a grid point is defined as 


where 



“j,k * ■ l yV' , n ) 


(j-l)AC 

C_j. 

II 

H 

2 , . . . , JM 

(k-l)An 

k = 1, 

2, . . . , KM 

n At 

n = 0 , 

1, . . . 


Referring to figure 3.1(b), we see that this definition puts 5 = 0 along the 
left edge and n = 0 along the bottom of the grid, and restricts the number of 
moving points to 2 ^ j ^ JH; 2 £ k £ KH. 

The transport equations and the equation for the stream function are quite 
different in character. Both the vorticity and energy equations are parabolic, 
time-dependent partial differential equations; whereas the Poisson equation for 
the stream function is elliptic and explicitly independent of time. This leads 
to the following generally accepted computational pattern. Given the values of 
all dependent variables over the mesh at the beginning of a time interval: 
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1. Advance the vorticity and energy one time increment (either implicitly 
or explicitly). 

2. Using the newly evaluated vorticity, compute the stream function by 
solving the Poisson equation. 

3. For increased accuracy, cycle step one and two in a conventional 
predictor-corrector process. 


VORTICITY AND ENERGY EQUATIONS 


Finite-Difference Scheme 

The finite-difference approximations to the transport equations are based 
on the predictor-corrector scheme developed by MacCormack (ref. l). The dif- 
ference formulations are in conservative form and are accurate to second order 
in both space and time, having an error proportional to a third space deriva- 
tive in dispersion and a fourth space derivative in dissipation (see ref. 2). 


Predictor for the vorticity: 
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(3.2a) 


Predictor for the energy: 
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Corrector for the vorticity: 
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(3.2c) 
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Corrector for the energy: 


■Ti + - ( ™W ] - It - ( ™W 


+ ? 7T7TT (H , . - 2H. . + H. , ) + - / T v_ (H. 2H. , + H. , )\ 

Pr(A?)2 j+1 ,k j ,k j-1 ,k Pr(Art)2 j ,k+l j 5 k j ,k-l \ 


(3. 2d) 

In these expressions the convective terms are differenced forward and backward 
in space for the predictor and corrector equations, respectively. In practice, 
however, this procedure is programmed to be cyclic so that all four possible 
combinations of forward-backward differencing are used in four successive time 
steps. In the intermediate step U and V are obtained via the solution of the 
Poisson equation for ¥ from ft. Similarly, in the final step U n+1 and V n+1 
are obtained by the solution of the Poisson equation for from a) n+1 . 

Notice that the diffusion terms are central differenced in both the predictor 
and corrector. The term Aft in equations (3.2a) and (3.2c) is evaluated 
according to Aft = (dft/dx)Ax where dft/dx is a specified time-dependent 
function. 


Stability 

The stability of a general set of nonlinear difference equations is 
usually estimated by studying the stability of an "equivalent" set of linear 
difference equations. We accept this philosophy and point out further that not 
only can the stability of the linear difference equations be established but 
their complete analytic solution can be determined in a straightforward manner. 

Consider, for example, the simple one -dimensional diffusion equation 


M. 

dT 


y 


8 2 cj) 


A well-known difference version is 


2 - 2*“ + .) 
z . 1+1 .1 . 1-1 


,n+l ,n . At 

b h V (A?)2 -J + l ' j • j- 


which can be expressed in vector, matrix notation by 

d> = Bd> + f 
x n+i 

or 

(EI - B)$ n = f 


(3.3) 


(3.4) 


(3.5a) 


(3.5b) 
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where E is the displacement operator, and f contains the boundary condi- 
tions. Note that when there is no space index, the time index is written as a 
subscript rather than a superscript. This is done in this section to clarify 
exponentiation. If f is a vector of constants 1 the solution of equation 
(3.5) can be written (I being the unitary matrix) 

<f> = l c (X ) n + [I - B] -1 f (3.6) 

n j j j 

The c . are constants determined by the initial conditions , and X^ axe the 
values 0 of the scalar E, which are the assumed distinct roots of the character- 
istic polynomial equation 


P(E) e det (El - B) = 0 


(3.7) 


The matrix B and the vector f in equation (3.5) depend on the choice 
of differencing scheme and the boundary conditions. If they represent equation 
(3.4) with boundary conditions given at the two ends 5=0 and 5 = 1, they 
become 


B = T(g, 1 - 2g , g) 
f = (g<j>i , 0 , . . . , 


°> S *JM> 


T 


where g = pAx/A5 2 * and T 


is a square tridiagonal matrix defined by 


(3.8) 


T(d_i,do,di) = 


d 0 d x 0 

<3-1 <lo <3l 

0 d - 1 d 0 


0 0 
0 0 
0 0 


0 0 0 
0 0 0 


The eigenvalues of 


T(d_ 1 ,do,d 1 ) are 


do d l 
d_i d 0 




= do + 2/dqdIY cos 



3 = 1 , 


M 


(3.9) 


( 3 . 10 ) 


where M is the number of moving points. Notice that if B is independent of 
E, the roots of det(EI - B) = 0 coincide with the eigenvalues of B. Using 


1 That is, if the boundary conditions are independent of x. When such is 

not the case the particular solution is more complicated, but can be found (see 

refs. 3 and U). 



2 fixed points I, JM 
JH-I moving points 2 < j < JH 



the normalized mesh spacing shown in sketch (a), we 
can find the eigenvalues of B in equation (3.8); 
and using equation (3.6), we can write the analytic 
solution to equation (3.4) as 


j= i 

t 




2 


3 


JH JM 

t 

£= l 


Sketch (a) 


JH-1 

<}> = ][ c.[l - 26 + 26 cos ( jir/JH) ] n + [I - B] _1 f 

~ n J=1 ~J 

(3.11a) 

under the imposed boundary conditions . Equation 
(3.11a) has the alternative expression 



1 - 


4y At .9 ,1tt A£ 
sin 2 


n 


+ [I - B ] - 1 f 


(3.11b) 


Clearly, as n increases the solution of a set of linear difference equa- 
tions can grow unboundedly if the absolute value of any Aj in equation (3-6) 
is greater than one. Two remarks should now be considered. First, it is pos- 
sible in an analytic construction to set to zero all the elements in the vector 
Cj multiplying a given A j . The magnitude of such A^ could be greater than 
one and a bounded solution would result as n + 00 . This behavior is basically 
related to a saddlepoint problem (ref. 5). Second, it is possible that all of 

the | A j | are less than one but they have a structure such that a computed 

solution will appear to be growing in an unbounded manner as the first several 
time steps are calculated. This behavior is discussed in reference 6, p. 152. 
Nevertheless, equations such as equation (3.1l) are the exact solutions of the 
linear difference approximations and they tell us precisely (except for round- 
off error) what a computer would calculate after any number of time steps. 


The stability of difference equations is generally not (for an exception 
see ref. T 5 p. 222) viewed in the above light. More often an amplification 
factor is developed along lines similar to those introduced by Von Neumann, as 
discussed in reference 8. Thus the term 

ik^C ik n 

4>(t,€,n) = <f>(x)e ^ e n (3.12) 


is introduced into the difference equation and the ratio A = <f>(x + At)/<J)(t) is 
determined. This ratio is referred to as the amplification factor and the 


values 
is , of 
finds 

or 


k^ and k^ are called the wave numbers 
course, that | A | <1. 


The condition for stability 
Applying this technique to equation (3.4), one 


X = 1 + B( 


ik ? A5 


- 2 + 


-ik AC 
e ? ) 


(3.13) 


A = 1 - 26 + 26 cos k- AC 


This appears to be the same as the term inside the parenthesis in equation 
(3.11a), but it is not because the arguments of the cosine terms are different. 
We can make the two approaches identical, however, if we reexamine the discus- 
sion under equation (3.7). 
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The solution given "by equation (3.1l) depended on the nature of the matrix 
B, which in turn depended on the boundary conditions. Suppose, instead of fix- 
ing the conditions at 5 = 0 and E, = l,*we require that the solution he peri- 
odic. In such a case equation (3.8) would become 


B = T (0, 1 - 20, 0) 

f = ( 0 , 0, . . . , 0, o) T 


where T 
~P 


is a square "periodic-tridiagonal" matrix defined by 


(3.14) 


Tp ( d_ i , do ,di ) 


d 0 dj 0 

d_i d 0 d : 

0 1 do 


0 0 0 

di 0 0 


0 d_! 

0 0 

0 0 

do di 

d- 1 do 


(3.15) 


the eigenvalues of which are 


Aj = do + ( d_ ! + di)cos + i ( d_ ! - d^sin - 1 


j— 1,2, . . • * M 

( 3 . 16 ) 


Using equations (3.1^) and (3.16), we see that the exact analytic solution to 
equation (3.^) is 


JH-1 

cf> = l c . 

J = 1 J 


1 - 23 + 23 cos 


2tt(,1 - 1) 
JH-1 


(3.17a) 


when the boundary conditions are periodic. This has the alternative form 


JH-1 

= l 

j=i 


~ c j 


i - 


4p At 
A52 


sin 2 (tt ( j - l)A?) 


(3.1Tb) 


^jh 1 Periodic boundary 
4*2 z J conditions 


-|‘ £ r 

• • • 

) = l 2 3 

t 

€ = 0 


= 


(JH-1) 


JH JM 

t 

e= I 


where A£ = l/(jH-l) for a periodic mesh (see 
sketch (b)). The term in the brackets in equation 
(3.17) corresponds to equation (3.13). 

The analytic solutions to linear difference 
equations with periodic boundary conditions are 
especially easy to find because the eigenvalues in 
their B matrix can be readily determined. In 
fact, the simplest way to find these eigenvalues is 


Sketch (b) 



often by the separation of variables technique employed in the discussion of 
equation (3.12). Consider, for example, the model equation that couples con- 
vection and diffusion 


9<f) 

3x 




(3.18) 


The finite-difference approximation equivalent to the one used on equation 
(3.1) to form equation (3.2) gives 


»j * + 


(3.19a) 


♦f* = I [ +j + } j * “<*J - ♦j-l ) + 6< Vl - 2 +J + (3 - 19b) 


where a = c At/A£, often referred to as the Courant number. Substituting the 
predictor equation (3.19a) into the corrector equation (3.19b) gives 


,n+l 




- f <c 2 - - o 


+ 4 (+“ - H n n + 6<d - H n + <i> n ) 

2 VV j+2 ? j + l \j-l V j-z' 


( 3 . 20 ) 


The B matrix for this difference equation with periodic boundary conditions 
is the "periodic-pentadiagonal" matrix 


d 0 

d l 

d 2 

0 

0 

. . 0 

d_ 2 

d_T 

<3-1 

d 0 

d l 

d 2 

0 

0 

0 

a -2 

<3-2 

<3-1 

d 0 

d l 

d 2 

0 

0 

0 


P (d_ 2 » d -l5 d 0 5 ( il 5 < i2) = 
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0 

0 
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0 

0 
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<30 
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d 2 

0 

<3-2 

<3-1 

d 0 


where 


(3.21a) 



d_ 2 = 


ctg/2 

+ 

B 2 /2 



II 

rH 

J 

a/2 + a 2 /2 + 

6 - ct3 

- 

2B 2 



d 0 = 1 

- a 2 - 

26 

+ 

3B 2 


> (3.21b) 

dl = - 

a/2 + a 2 /2 + 

6 + aB 

- 

2B 2 



II 

CM 

n3 


- aB/2 

+ 

6 2 /2 

J 



and its exact analytic solution is 


JH-1 

Sj'V" 


where 


x j 


= 1 + (a 2 + 23)< 


COS 


2tt(,i-1) 

JH-1 


- 1) + 2g 2 \cos 


2 TT ( .j — 1 ) 
JH-1 


- 1 


-i \ (a - 2a3)sin 


2tt(,1-1) 

JH-1 


+ aB sin 


4tt(,i-1) 

JH-1 


J = 1, 


(3.22a) 


JH-1 

(3.22b) 


these being the eigenvalues of equation (3.21a) with elements defined by equa- 
tion (3.21b). Although these eigenvalues are complex, each complex has a 

conjugate so the summation (3.22a) results in real numbers. 

It can be shown that | Aj | has extrema at j = 1 and at j = 1 + (jH-l)/2. 
Notice that | A i | =1 for all a and 3. Designate the value of Aj at 
j = 1 + (jH-l)/2 by A c , and one can write 

| Aj = |l - 2a 2 + UB(23 - 1) | (3.23) 

The stability boundary is formed by finding the values of a and 3 for which 
| A c | <1. This boundary is formed by two lines: one, where A c < 1 , which 

results in 

2g(2B - l) - <x 2 < 0 (3.24) 

and the other, where A c > -1 , which results in 

a 2 + 20(1 -20) < 1 (3.25) 


The a, 3 combinations that produce stability for equation (3.19) with 
periodic boundary conditions are shown in figure 3.2. The rectangle for which 
0 < | a | < /3/2 and 0 < 3 < 1/2 represents the practical stability boundary. 

Note that 3 = y At/A^ 2 is taken to be positive since only positive diffusion, 

^9 




viscosity 5 or thermal conductivity 
coefficients are considered. On the 
other hand a = c Ax/Ag can be posi- 
tive or negative depending on the 
direction of the wave velocity. Note 
that the effect of coupling diffusion 
to convection is to lower the maximum 
generally allowable value of the 
Courant number from 1 to 0.866. It 
follows from the above that the stabil- 
ity of the method represented by equa- 
tion (3.19) when applied to periodic 
boundary conditions requires that At 
be chosen so that 


Figure 3.2.- Stability boundaries for equations (3.19) 
with periodic boundary conditions. 


At < 


both 


Mi 

2y 


and 


/3 _AC 
2 I c I 


(3.26) 


In the more important two-dimensional case the concepts described above 
remain the same but the algebra becomes more involved. The partial differen- 
tial equation coupling diffusion and convection for two-dimensional flow is 


11 = _u _ 

3t 35 





(3.27) 


where U, V, and y are constants if the equation is to be linear. If the 
predictor-corrector sequence, defined for one- dimensional flow by equation 
(3.2), is applied to equation (3.27), and periodic boundary conditions are 
imposed, detailed calculations show that the stability boundaries for | X c | 
are the same as those given by equations (3.21+) and (3.25) except that 


•- Nfr* MS 


e = y At 


1 + 1 


A? 2 An 2 




( 3 . 28 ) 


This means that the difference equations with periodic boundary conditions are 
stable if 


Ax £ 


both 


_1 ( A5 2 An 2 \ 
2y 1 A5 2 + An 2 J 


and 


/3 A5 An 

2 l|u| An + I V I AC 


(3.29) 


In the application to physical problems the effectiveness of the 
predictor-corrector method, described by equation (3.19) for the linear and 
equation (3.2) for the nonlinear cases, respectively, depends on the relative 
magnitudes of the diffusion coefficients (i.e., the viscosity and thermal con- 
ductivity) and some representative or average value of the velocities. If 
there is a high degree of diffusion (y large) then the first term in equation 
(3.29) would be the smallest and the time step size would be forced, by reasons 
of stability, to be unnecessarily small for a given accuracy. In such cases an 
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implicit differencing of the diffusion terms is suggested (refs. 9, 10). If, 
on the other hand, the diffusion is low relative to the magnitude of the (ave- 
rage) velocity, the physical requirements associated with the domain of depen- 
dence limit the size of the time step, and the explicit method described is 
entirely adequate. 

In the particular case regarding the convective mixing of oxygen at high 
pressure and low temperature, both the viscosity and the thermal conductivity 
are so low (ch. 6) that the explicit method is entirely satisfactory from the 
viewpoint of stability. 


Accuracy 

The accuracy of the differencing scheme given by equation (3.19) can be 
estimated in several ways. By accuracy we mean, of course, how well it repre- 
sents the basic partial differential equation (3.18). We choose here to dis- 
cuss this accuracy by deriving the "modified partial differential" equation, 
that is, the differential equation actually represented by the difference equa- 
tion (see ref. 2) and comparing this to the basic partial differential equa- 
tion. One can derive the modified partial differential equation for equation 
(3.19) by expanding each term in a Taylor series about the point represented by 
(n,j). Time derivatives of higher than first order are eliminated by repeti- 
tive use of equation (3. 20) itself. The result of such a procedure gives 


= ~ C(t> E, + - - c 2 (At) 2 ]<)>^ 

-^q^[(A 5) 2 - c 2 ( At ) 2 ] - ^[(AS) 2 - 6c 2 (At) 2 ]|^ cc? + . . . (3.30) 

The coefficient of the third derivative represents, to the lowest order, the 
dispersion or phase error, and the coefficient of the fourth derivative repre- 
sents, again to the lowest order, the dissipation or diffusion error. 

Note that the error in dispersion is due only to the finite differencing 
of the convection term (given by the coefficient of c in eq. ( 3 . 18 )) and is 
of order (A£) 2 , (At) 2 . The error in dissipation or diffusion is caused by the 
finite differencing of both the convection and diffusion terms, orders being 
(At)(A£) 2 , (At) 3 , and (A£) 2 , (At) 2 , respectively. When the magnitudes of the 
viscosity and thermal conductivity coefficients are small, as they are in our 
intended physical application, the principal error is in dispersion and is pro- 
portional to (c/6)[(a?) 2 - c 2 (At) 2 ]. In the practical problem, c represents 
a local velocity and the magnitude of the velocity varies throughout the tank. 
One can show that the term (c/6)[(A^) 2 - c 2 (At) 2 ] is zero when c At/A£ is 
either zero or one, and a maximum in that interval when c At/A£ = 1//3. Since 
|c At / A£ | max < 1 for stability, the maximum error in dispersion given by the 
lowest order error term in the modified partial differential equation (3.30) is 
-(c A£ 2 /9)<(>£££ if numerical stability is the only bound on the time step size. 
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To get less dispersive error with the differencing scheme given in equation 
(3.19) one must choose a time step such that c At/A? < 1/^3. 


POISSON EQUATION 


The finite-difference approximation to the Poisson equation (3.1c) for the 
stream function may he written with second-order accuracy hy the well-known 
five-point formula as 


, - 2¥. . + 

j ? k - . 1 14 

(An ) 2 

j = 2, . 


AA±L + 


- 2 ¥ 


JM-l 


k = 


■i ,k 
(A?) 2 

2 , . , 


+ ¥ 


Jil _ 


KM-l 


—to . , 
J 


(3.31) 


Given a (JM-2) x (KM-2) internal distribution of vorticity and 2(jM + KM - 2) 
boundary values for the stream function as shown in figure 3.1(b), we wish to 
find the (JM-2) x (KM-2) interior values of ¥ that satisfy equation ( 3 * 31 ) . 
We can find such a solution by iterative techniques such as SOR and ADI (ref. 
ll) , but since the equation is linear we can also solve it by direct inversion 
of a set of linear, simultaneous algebraic equations. These so-called direct 
methods require far less computing time than the iterative ones when used to 
find solutions to equation (3.3l). The particular one used here was developed 
by Buneman (ref. 12) and is referred to as the double cyclic reduction method. 
It is an extension of the odd/even reduction scheme (refs. 13, lU). 


Odd/Even Reduction of a Tridiagonal Matrix 


Consider the tridiagonal set of algebraic equations 


a -1 0 

-1 a -1 

0 -1 a 

0 0-1 
0 0 0 

0 0 0 

0 0 0 


0 0 

0 0 

-1 0 

a -1 

-1 a 

0 -1 

0 0 


0 

0 

0 

0 

-1 

a 

-1 



(3.32) 


Multiplying the odd rows by a and adding the adjacent rows to it gives 


— — 


— — 



a 2 -2 -1 0 


u 3 


af 3 + f 2 + f4 

-1 a 2 -2 -1 


u 5 

= 

af 5 + fit + fe 

0 -1 a 2 -2 

L J 


_ U7 -J 


7 + f6 + 


(3.33) 
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which is obviously recursive if the number of rows equals 2 -1 where L is an 
integer. If we set 

1 ^ = a , a^ +1 ^ = [a^] 2 - 2 (3.34a) 

fO) _ f f (^+l) _ ( 3 . 34"b ) 

n n n n „£-l 


n-2 


n+2 


-L-l 


and let M represent the midpoint in the U vector (M = 1 + 2 ) , we have 

after L-l recursions 


&(l)u m = 4 L) 


(3.35) 


This is solved for u^ and a backward recursion gives the completed solution. 

One of the advantages of this scheme is the efficient use of core storage. 
Only one array is needed in the calculation, it is initially filled with the 
elements of f , which are then overwritten and replaced by the values of u in 
the solution. Table 3.1 shows how an array 2 of IT is overwritten by the for- 
ward recursion, at the end of which every other value of f has been overwrit- 
ten. Table 3.2 shows the backward recursion. The final value of ug is 
computed by equation (3.35), then u 5 and u 13 are found, and so on, until all 
the u from 2 to 1 6 are evaluated and have replaced the original f 2 , . . . , 

f 16 - 


The same concept can be applied to find the direct solution of the two- 
dimensional equation (3.3l). When equation (3.3l) is multiplied by (A£) 2 and 
written in matrix form, we obtain 


A -I 

-I A -I 

-I A -I 

0 


0 




r n 




CNJ 




*3 


£3 


l 




-JM-1 


(3.36) 


where A is a (KM-2) x (KM-2) tridiagonal matrix with -2[l + (A$/An) 2 ] on the 
diagonal and (A?/An) 2 on each side; I is the (KM-2) x (KM-2) identity matrix; 

is the vector of stream functions in the jth column of the grid; and f-j 
is the jth column of vorticities and boundary values. Note that JM must be 

2^+1 and KM must be 2^ 2 +l where Iq and L 2 are integers. 


2 

L = h giving 2 4 - 1 = 15 moving points, but a word at each end is 
reserved and set to zero to simplify the backward recursion; see table 3.2. 
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TABLE 3.1.- FORWARD RECURSION IN ODD /EVEN REDUCTION. COLUMN 
ON LEFT IS ORIGINALLY STORED, THEN SUCCESSIVELY USED AND 
OVERWRITTEN BY COLUMNS 2, 3, AND k 



TABLE 3.2.- BACKWARD RECURSION IN ODD /EVEN REDUCTION. COLUMN ON 
LEFT IS ORIGINALLY STORED, THEN OVERWRITTEN BY COLUMNS 2, 3, k, 
AND 5 


fl ^-0.0 

f 2 

f 3 

f 4 

f 5 

*6 

f 7 

fg ^ f 9 /a( 4 ' 

fio 
fll 
f 12 
f 13 
f 14 
f 15 
fl6 

f 1 7 0.0 


(f5+fl+fg)/a( 3 ) 


(fl 3 +f 9 +fl 7 /a^ 3 ) 


(f 3 +f!+f 5 )/a( 2 ) 

( f 7 +f 5 +f 9 ) /a^ 2 ^ 

( f 1 i+fg+f i 3 ) /a^ 2 ) 

( fl 5 +f 1 3 +f, l 7) / 


(f 2 +f 1 +f 3 )/a(^ 

(f4+f 3 +f 5 )/a( 1 ^ 

( f 6 +f 5+f 7) / a( 1 ) 

(f 8 +f 7 +f 9 )/a( ^ 

( f 1 o + ^ 9 + ^ 1 1 ) / a ^ 1 ^ 
(fl2+fn+fi 3 )/a( ^ 
(fl4+fl 3 +fi 5 )/a^ 1 ^ 
( f 16 +f 15 +f 17)/ a ^ ^ 






The process of odd/even reduction described for a simple matrix such as 
that given by equation (3.32) can now be applied to a block matrix such as that 
given by equation (3.36). The concept is simplified if we refer at once to 
tables 3.1 and 3.2 and recall that matrix and algebraic arithmetic are identi- 
cal in addition and differ in multiplication only on commutivity. Consider, 
first, table 3 . 1 . Wote that all operations are formed by simple additions and 
multiplications of the data existing in the single f array at any given stage 
in the calculations. The same is true in the two-dimensional case, except that 
the operations are on the column vectors f. For example, for the two- 
dimensional case, in table 3.1 the first entry of the second column would be 
replaced by A .( 1 )f 3 + f 2 + fi* , which is a tri diagonal multiplication into a 
vector and two vector additions . 


The two-dimensional equivalent of equation (3.3^a) is formed by the matrix 
definitions 


( 1 ) _ 


u+i) _ A U) fl U) 


= A , A' ' = A A - 21 

The essential part of this relationship is that A^ +1 ^ 
to a product of tridiagonal matrices. Thus 


(3.3T) 

can always be reduced 


A.^ 1 ^ = A 

A^ = (A - v/2 i) (A + /2I) 

a^ = (a - J2 - J2 . i)(a - Jmi i ) (a + /2 - 1) (a + v4T+~7i i) 


A^ = (A - X^l) (A - X 2 l) ... (A - X^X ) 

,£-l 


\ (3.38) 


k = 2 


X. = 2 cos 
<3 


( 2,j -1 ) 7T 

2 k 


j — 1 j 2 , . 


, k 


(£) 

In general, A is the product of 2 tri diagonal matrices that have constant 
entries along the two off diagonals and differ only by the constant Aj along 
the diagonal. If we write a subroutine that will perform such a simple opera- 
tion, we can apply it twice in succession to form each of A' 2 'f 5 5 A^'fg, and 
( 2 ) ( 3 ) 

A f 1 3 , and four times to form A fg in the two-dimensional equivalent of 

table 3.1. 


The remarkable property of the process being described comes in the two- 
dimensional equivalent of the procedure outlined in table 3-2. The matrix 
equivalent of an algebraic division is a product using the matrix inverse. 

This would greatly complicate the extension of the backward recursion if the 
factorization shown in equation (3.38) did not exist. Since such a factoriza- 
tion does exist, however, the actual computation of u^ in the two-dimensional 
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equivalent to equation (3.35) 


^ (I \ ■ 4 L) (3 - 39) 

is found by performing 2 L 1 tri diagonal inversions (each identical to the 
scalar manipulations shown in tables 3.1 and 3.2 on the data in f^. For exam- 
ple, the two- dimensional equivalent of (f 3 +f i+f s)/a( 2 ) in table 3.2 is 


(A - /2 I)" 1 (A + v^2 I )” 1 ( f 3 + fi + f 5 ) (3. HO) 

and is computed by two successive calls to a tridiagonal inversion subroutine, 
each call having a different constant for the diagonal entry but both with a 
fixed constant for the two off-diagonal entries. 

This describes the concept of one- and two-dimensional odd/even reduction 
of the Poisson equation. In application to the two-dimensional case (eq. 
(3.3l)) the matrix products formed by the forward recursion can lead to very 
large floating point numbers. A clever way to avoid this numerical complica- 
tion has been devised by Buneman and a FORTRAN program that computes the 
Poisson equation for Dirichlet boundary conditions using Buneman T s version of 
odd/even reduction is listed in reference 12. 


BOUNDARY CONDITIONS 


The boundary condition for the energy equation (3.1b) is satisfied by 
specifying either the temperature or the heat flux at the wall. The boundary 
condition for the Poisson equation (3.1c) is satisfied by setting ¥ = 0 along 
the edges. This corresponds to the condition of no mass flow through any por- 
tion of the wall. The boundary condition for the vorticity equation (3.1b) is 
somewhat more subtle and is described below. 

For a viscous fluid, the additional condition of no slip at the walls is 
met if the normal stream- function derivative (or tangential velocity) vanishes 
at the walls. This boundary condition is satisfied in the vorticity equation 
(3.1a) in the following way. First, notice from equation (3.1c) that along a 
horizontal wall 



0) = -Y 

(3.1+la) 

and along a vertical wall 

nn 


“ = 

(3.1+11) 


since ¥ itself Is zero along the walls. Next consider a Taylor 1 s series 
expansion for the stream function at a point on the wall. Using it to find the 
value of ¥ at an interior point next to the wall, one finds (taking the bot- 
tom wall as an example) 
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(3.42) 
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Again 'Fj- 1 must be zero since it is the value of ¥ at the wall, and now we 
can apply’the no-slip condition by also setting the tangential velocity at the 
wall (3f/9n)^ 1 equal to zero. This reduces to a boundary condition for o>; 
namely, at the wall a> is given by 
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3,1 



(3.43) 


where ^ as ^ een calculated in a previous predictor or corrector sequence. 

Similar formulas are also derived for the other walls. Formulas of higher 
order may be obtained by expanding about deeper interior points; for example, 
the second-order approximation for points ( j ,l) is found by simultaneously 
solving the expansions 


AH' . o = A 
J >2 


j »i Un 


. < 4 ” )+ ! pj . < 4 ”) 2 + 


2 + i 
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6 l 8n 3 


( An ) 3 + 0(An) 4 
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(3.U5) 


Singularities exist in the analytic solutions for the flow at the corners 
of the box, and it is possible for a corner to lie outside the radius of con- 
vergence of a Taylor series expansion about one of the neighboring points. In 
fact, the calculations did show instabilities near the corner points when the 
higher order equations were used, and it is believed they developed for the 
reason just described. As a result, the lower order equations, typified by 
equation (3.43), were used for the boundary condition on the vorticity through- 
out the calculations. The sacrifice in accuracy is not considered important 
for the nature of the mixing problem involved. 


COMPUTATIONAL PROCEDURE 


To begin the convection calculations, the initial velocities are set to 
zero and the temperature distribution specified. The subsequent computational 
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procedure is described by the following steps: 

1. Calculate the predicted vorticity and temperature distributions from 
equations (3.2a) and (3.2b). 

2. Calculate the corresponding stream function by Buneman's double cyclic 
reduction method. 

3. Apply the wall boundary conditions for the vorticity from equation 

( 3 . 43 ). 

4 . Calculate the corrected vorticity and temperature distributions from 
equations (3.2c) and (3. 2d). 

5. Calculate the corresponding stream function as in step 2. 

6. Apply the wall boundary condition as in step 3. 

The calculations for one time step are now complete. The procedure continues 
until the desired number of time steps are complete. 


CONCLUDING REMARKS 


A finite-difference predictor-corrector scheme has been described for cal- 
culating convection in a rotating tank. The assumption is made that the sta- 
bility bounds and the order of accuracy of a set of nonlinear difference equa- 
tions can be estimated by studying an "equivalent" set of linear model 
equations. Such being the case, it is pointed out that not only can the sta- 
bility of the linear difference equations representing convection be estab- 
lished, but the complete analytic solution of them can be written for the 
n(th) time step without ever going to a computer. Furthermore, such an ana- 
lytic solution to a set of linear difference equations with periodic boundary 
conditions predicts stability boundaries identical to those obtained from an 
amplification factor analysis. The accuracy study for these convective equa- 
tions is based on a study of the modified partial differential equations, and 
it shows that the approximations are characterized by second-order errors in 
both dispersion and dissipation. 

Finally, a fast, efficient, and direct numerical solution to Poisson ! s 
equation is described. It is based on double cyclic reduction techniques 
introduced by Buzbee, Golub, and Nielson in reference 13. 
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4 . CALCULATION OF THERMODYNAMIC PROPERTIES OF OXYGEN 
NEAR THE CRITICAL POINT 

Walter A. Reinhardt 





SUMMARY 


Methods are developed for the versatile and efficient evaluation of the 
thermodynamic properties of cryogenic oxygen. The semi empirical equations of 
state of Stewart are used in this discussion. An arbitrary choice of indepen- 
dent thermodynamic variables within a limited set is allowed. Comments on pro- 
cedures that lead to an expansion of the variable set are given. A method for 
the efficient rapid evaluation of integrals representing volume averages of 
spatially varying thermodynamic quantities is described. A number of graphs 
are given that show state relations plotted in different ways to demonstrate 
the methods used. In addition to being instructive, these curves serve as a 
valuable ready reference on cryogenic oxygen properties. 


INTRODUCTION 


The principal objective of this work is to formulate the thermodvnamic 
properties of cryogenic oxygen to permit their efficient evaluation by elec- 
tronic computers. The resulting computer programs are used in conjunction with 
other codes in the simulation study of the behavior of gases stored in tanks of 
maneuvering space vehicles. Comments on the complementary numerical codes 
along with discussions on the complete problem and background surveys related 
to this report are given in chapter 1 of this report. It is worthwhile to 
orient the reader, however, on the material presented in this article by 
briefly reviewing a few pertinent considerations that necessitated the availa- 
bility of accurate thermodynamic properties that can be rapidly computed. 

In the case of the Apollo spacecraft, the onboard cryogenic storage system 
provides gaseous hydrogen and oxygen (oxygen is of particular interest in this 
study) used in the electrical power system and environmental control system of 
the command and service module. Fluids are stored in the system at tempera- 
tures and pressures slightly greater than critical where the density is so 
large that the fluid is not distinguishable as being either liquid or gas. 
Cryogenic storage allows spacecraft minimum volume and weight requirements to 
be met and ensures single-phase delivery of the fluids (ref. l). To maintain 
a constant pressure within the tanks, as required by the electrical power and 
environmental control systems as oxygen is drawn from the system, heat is added 
by means of contained electrical heaters. The added heat can remain in the 
neighborhood of the heaters, however, in the absence of conductive or convec- 
tive mixing; there then result large thermal gradients ("thermal stratifica- 
tion") within the stored fluids. Because oxygen has a very low thermal conduc- 
tivity, conduction currents are negligible. Also, convective mixing processes 
do not take place in the absence of any spacecraft accelerations in the other- 
wise gravity- free environment. If particularly severe thermal stratification 
occurs, a spacecraft maneuver (such as a midcourse correction on a flight to 
the moon, for example) could lead to mixing that would then result in a large 
abrupt pressure decrease or "pressure collapse" (ch. 1 and ref. 2) with atten- 
dant possibility of catastrophic equipment failure. Prior to Apollo ik , fans 
were included within the storage tanks to ensure that the fluid was mixed and 
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uniform temperatures maintained. The subject study of this report is that of 
investigating the use of vehicle maneuvers (in particular , rotations) as a 
means toward the maintenance of uniform temperatures and pressures. The nomi- 
nal operating pressure requirements of the cryogenic gas storage system are 
900 ±35 psia (62 ±2,4 atm) or, stated in a slightly different manner, that the 
pressures must not deviate by more than about 4 percent from the mean stated 
value. This suggests that accurate relations are required for gas properties 
if reliance is to be placed on interpretations gained from numerical computa- 
tions, The semiempirical thermodynamic properties published for cryogenic oxy- 
gen by Stewart (ref. 3) and by Weber (ref. 4) satisfy the accuracy requirements . 

One complication that occurs in the analytical description of fluid prop- 
erties near the critical point is that the properties require two independent 
variables with both variables showing strong nonlinear behavior. In this 
respect, the relatively simple descriptions such as provided by van der Waals r 
or Beattie-Bridgeman’ s equations of state (see, e.g., ref. 5 or 6) are not at 
all satisfactory. In fact, the required state equations can be quite compli- 
cated, as shown in both Stewart’s and Weber’s papers (refs. 3 and 4). An addi- 
tional problem that occurs with such relations is that should a different 
choice of independent variables be required, considerable computer time can be 
expended in the application of conventional inversion methods. 

In the discussion that follows, particular emphasis is placed on formulat- 
ing the thermodynamic properties of cryogenic oxygen in such a manner that 
flexible and time-wise efficient evaluations are possible. By flexible it is 
meant that an arbitrary choice of independent variables is allowed. A minimiz- 
ation of redundant arithmetic operations leads to an efficient computer program. 

It was found that a "modified virial" representation of Stewart’s equa- 
tions leads to the least computational redundancy. Finding zeros by means of 
successive linear interpolations (ref. 7) provides the means for obtaining 
rapid function inversions. The thermal equation of state and the caloric equa- 
tion of state are plotted in several different ways to demonstrate the flexi- 
bility of the procedures described. 

A method is also described that is particularly valuable for the rapid 
evaluation of volume integrals of thermodynamic spatially varying quantities 
(e.g., temperature) where the integrals represent volume averages. An impor- 
tant concept of this method is that accuracy can be traded to obtain computing 
speed. Thus, by sacrificing the numerical accuracy beyond that specified by 
the semiempirical formula, one obtains significant savings in computational 
time. The method is particularly well suited for certain classes of numerical 
simulation studies, since evaluation time is nearly independent of the number 
of spatial mesh points . 

The author wishes to acknowledge several helpful conversations with Prof. 
Richard B. Stewart, Mechanical Engineering Department, University of Idaho, 
Moscow, Idaho, and with Dr. Loyd A. Weber, Cryogenics Division, National 
Bureau of Standards Laboratory, Boulder, Colorado. 
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THERMODYNAMIC QUANTITIES 


There are several accurate representations for the thermal equation of 
state for cryogenic oxygen. In one case (ref. U) 5 accurate experimental 
results in tabular form are published that encompass a broad range of values of 
pressure and temperature (5^-6° T < 300° K; 1 < p < 330 atm). A value of a 
state variable , say pressure, can be made available on an electronic computer 
by the use of a table-lookup and interpolation procedure once the other pair of 
variables , density and temperature, are specified. In the second case (ref. 3), 
the state variable is represented by a complicated analytical function of two 
independent variables. Specification of values for the two independent vari- 
ables and evaluation of the function are then necessary to obtain a value for 
the state variable. 

A purpose of this work is to design a flexible procedure for obtaining 
thermodynamic state variables - that is, a procedure that allows the rapid 
determination of an arbitrary variable given any pair of other variables. To 
expedite this task, it was decided to minimize the coding problems encountered 
by representing the state equation as a complicated analytical relation. As a 
result of interpolation problems , table-lookup procedures that depend on two 
independent variables often are inapplicable to general usage. For example , 
when derivatives of thermodynamic quantities are required, numerical differen- 
tiation of table quantities can result in an unacceptable lack of smoothness . 
Although one can develop procedures to overcome these difficulties, they are 
avoided by the use of the analytical relation. 

In the discussion that follows the analytical representation given by 
Stewart in reference 3 is used. It will also be shown that by modifying 
Stewart T s original formulation to one similar to the virial state equation 
(see, e.g., ref. 6) one can significantly reduce computer evaluation time. 

One should note that although Stewart’s formulation was used, much of the 
content of this chapter is of a general nature and therefore applies indepen- 
dently of this particular choice. 


Quantities Defined Explicitly 

This section concerns the various thermodynamic quantities that are 
defined explicitly in terms of representations containing the independent vari- 
ables. Particular emphasis is given to formulations of the state equations 
where redundant arithmetic operations are held to a minimum; hence, evaluation 
time is also minimal. 

The particular state relations to be discussed are the analytical formula- 
tion of the thermal equation of state p = p a (T,p), the equivalent "physical" 
representation p = Pp(T,p), and the caloric equation of state e = U(T,p). 

Thermal equation of state (analytical representation ) — The equation of 
state that results from Stewart’s work is given by 
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p = p a (T,p) 


where 


RTp 


n 3 . n 4 . n 5 ) _ 2 


+ ^ ni T + n 2 + 52 + + T^/P 

+ ^n 6 T 2 + n 7 T + n 8 + + -^fjp 3 
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+n 24 exp [ n 2 7 ( T-T ) 2 ] p 28 v exp(n 2 6 'v 2 ) 


n 28 n 28 
V = p - p„ 


(4. la) 
(4.1b) 


In equations (4.l) p, T, and p are the variables representing pressure, tem- 
perature, and density, respectively. The symbol p a (T,p) denotes the right- 
hand side of the equation; the subscript c signifies critical point value. 

The quantity R is the universal gas constant. Values for the coefficients in 
the above equation, nj through ri 2 Q, are listed in the first column of table 4.1. 
Constants as well as some convenient conversion factors are listed in table 4.2. 
In reference 3, the coefficients were determined by the use of a least- 

squares procedure applied to experimental data provided largely by Dr. L. A. 
Weber of the National Bureau of Standards Laboratory, Boulder, Colorado. The 
measurements have since been improved by Weber (ref. 4). Differences between 
the latest measurements and the quantities computed by the use of equation 
(4.1a) are not significant (of order 0.1 percent) except in the neighborhood of 
the critical point where the density deviates by about 2 percent. The error 
for pressures and temperatures is less in this region. 

Several choices for the n^ given in equations (4.l) are listed in table 
4.1; they differ only in units. In the first two columns of the table, quanti- 
ties n^ are given both with and without a bar; the n± with bar are the 
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TABLE 4.1.- COEFFICIENTS FOR STEWART f S THERMAL EQUATION OF STATE 

(REF. 3) IN DIFFERENT UNITS 


i 

n i 

(p,atm; p ,g-mole/liter ; T,°K) 

n i 

(p ,g-mole/cc; T,°K) 

Cn i ] 

* 

( dimensionless ) 

1 

3.38759078D-03 

4.12851466D 01 

P" 1 

5.50454716D-01 

2 

-1.31606223D 00 

-1.60390749D 04 

T/p 

-1.381720910 00 

3 

-7.38828523D 03 

-9. 004229230 07 

T 3 /p 

-3.23827743D-01 

4 

1.92049067D 07 

2.34053474D 11 

t 5 /p 

3.51405940D-02 

5 

-2.90260005D 10 

-3.537^48190 l4 

t 7 / p 

-2.21722921D-03 

6 

-5.70101162D-08 

-6.94792010D-01 

(Tp 2 ) -1 

-1.91l60034D-02 

7 

7.968223750-05 

9.7H01019D 02 

P 2 

1.726314590-01 

8 

6 . 07022502D-03 

7.397886770 04 

T/p 2 

8.49721211D-02 

9 

-2.71019658D 00 

-3.30296280D 07 

(T/p) 2 

-2.45123930D-01 

.10 

-3.59419602D 01 

-4. 38030799D 08 

(T/p) 3 

-2.10038991D-02 

11 

l . 02209557D-06 

1.24564530D 04 

P- 3 

2.9524l697D-02 

12 

1.90454505D-04 

2. 32110 15 4d 06 

T/p 3 

3.55459957D-02 

13 

1.21708394D-05 

1.48328096D 08 

T/p 4 

3.02864054D-02 

l4 

2.44255945D-03 

2.97678886D 10 

(T/p 2 ) 2 

3.92722326D-02 

15 

1.73655508D 02 

2. 116369 30D 09 

T 3 /p 2 

1.014814650-01 

16 

3.01752841D 05 

3.67751334D 12 

(t 2 / p ) 2 

1.13936476D 00 

17 : 

-3.49528517D 07 

-4.25976365D l4 

t 5 / p 2 

-8.52721624D-01 

18 

8.86724004D-01 

I.08066567D 13 

T 3 /p 4 

9.21175040D-02 

19 

-2.67817667D 02 

-3.26393959D 15 

(T/p) 4 

-1.797654260-01 

20 

1.05670904D 05 

1.28782933D 18 

T 5 /p 4 

4.58284972D-01 

21 

5.63771075D-03 

6.87077425D 16 

(T/p 2 ) 3 

1.04114692D-01 

22 

-1.12012813D 00 

-I.36511926D 19 

T 4 /p 6 

-1.33656530D-01 

23 

1.46829491D 02 

1.789436050 21 

T 5/ p 6 

1.13200675D-01 

24 

9.9 8868924D-04 

3.057815290 03 

T/p 2n 28 

8.32903264D-03 

25 

-5.60000000D-03 

-5.6OOOOOOOD 03 

P -2 

-9.95505258D-01 

26 

-1.57000000D-01 

-3.94366040D 04 

p — 2n 2s 

-1.66253033D 01 

27 

-3.50000000D-01 

-3.50000000D-01 

T 2 

-1.46114912D-05 

28 

9.00000000D-01 

9.00000000D-01 

1 

9.00000000D-01 


Note: The last four digits beginning with D are representations for the 

exponent; for example, 9 . 00000000D-01 can also be written 9*10“* . 


TABLE 4.2.- CONSTANTS AND CONVERSION FACTORS 


R (Universal Gas Constant) = 0.0820535 liter- atm/ g-mole °K 

Critical point values for molecular oxygen (ref. 3): 
p c = 50. 14 atm 

T c = 154.77° K 

p c = 13.333 g-mole /liter 

e c = 742.2 J/g-mole 

Conversion factors: 

1 g-mole 0 2 = 31.9988 g (C 12 = 12.000 scale) 

1 liter-atm = 101.3278 J’(abs.) = 24.2179 cal = 0.0960417 Btu 
1 atm = 1.01325 *10 5 N/m 2 = 1.01325*10® dyne/cm 2 = 14.696006 Ib(wt)/sq in. 
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values listed by Stewart in his thesis (ref. 3) , and those without a bar. are to 
be used with the equations to be introduced here. These quantities (except 
ri 2 s through ^g) differ by a factor that is the universal gas constant. 

Stewart f s units for pressure, density, and temperature were atmospheres, 
g-mole/liter , and degree Kelvin, respectively (in Weber’s paper, ref. 4, the 
corresponding choice is MN/m 2 , g-mole/em 3 , and degree Kelvin). The coeffi- 
cients without a bar have units (n^) simply expressed as the ratio of powers of 
density and temperature (in Stewart’s units) as indicated by the symbols p 
and T listed on the right of the unbarred quantities in table 4.1. The con- 
version factors required for transformation of these constants are listed in 
table 4.2. In the last column of table 4.1 the coefficients, denoted as n*, 
are made dimensionless by appropriate ratios of the critical point values T c 
and p c . The dimensionless thermal equation of state that results can be use- 
ful also for the approximate evaluation of state properties of other nonpolar 
molecules (such as IT 2 * A, Xe, H 2 ) by introducing the tf prineiple of correspond- 
ing states" (ref. 6, p. 235). 

The analytical formulation, equation (4.1), is fairly complicated and is 
not in the most efficient format for frequent evaluation by an electronic com- 
puter. More useful is the modified virial formulation 1 given by 



=1+1 1 } ( T ) gf 1 } ( p ) (4.2) 

i=l 

Most conventional state relations can be cast in this form (see, e.g. , ch. 3 
and 4, ref. 6, and ref. 8). Although it can be argued that equation (4.2) rep- 
resents only a trivial reformulation of equation (4.1), it will become evident 
that worthwhile advantages are realized by capitalizing on certain features of 
the modified virial formulation. The quantities A^( x T(t) and g^( 1 )(p) are 
defined in the first columns of tables 4.3 and 4.4. The virial coefficients, 
the coefficients resulting from a power series expansion in density of the com- 
pressibility factor Z defined above, are readily found (ref. 6, p. 131, eq. 
(3.0-1) ) and are given by the following linear combinations of the A.C 1 ) 
quantities . 

B(T) = A^ 

C(T) = A ( 2 l) + A^ 1 ) 

D(T) = A^ 

E(T) = A^ 1 ) + A^ 1 ^ 


^ust and Gosman (ref. 9, p. 23l) have discussed the criterion required of 
a real-fluid thermal equation of state to ensure nonsingular state behavior as 
the ideal-gas limit of very low densities is reached. The thermal state equa- 
tions that can be placed in the modified virial form always satisfy their 
stated criterion provided that the functions g^( n )(p) are bounded as the den- 
sity becomes infinitesimally small. 




(4.3) 
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When the state equation is expressed in the modified virial form, approxima- 
tions are readily introduced that can lead to simplifications (e.g., for suf- 
ficiently small densities, say p < /| 10 5 /n 25 1 , only the lowest order virial 
term need be considered) . 

A more important aspect of the formulation is that the variable dependence 
is separated. That is, the temperature-dependent quantities are evaluated sep- 
arately from those that are density dependent. Use of this feature saves con- 
siderable computational time. For example, in the case where pressure and 
temperature are given and the density (implicitly defined) is to be found, as 
discussed in a later section, the required computer time is greatly abbreviated 
if the temperature quantities are evaluated only once. Further, since the 

coefficients A-P^T) are polynomials in temperature, it turns out that by 
separately evaluating and storing the individual terms (in parenthesis in 
table 4.3), one can use the same quantities during the subsequent evaluation of 
other thermodynamic quantities. For a given density and temperature, once the 
thermal equation of state is evaluated relatively little additional computer 
time is required to evaluate the other thermodynamic properties obtained by 
differentiation or integration of the thermal equation of state (internal 
energy, enthalpy, specific heat, and entropy). The significance of these com- 
ments will be particularly evident in the discussion that follows. 

Thermal equation of state (physical representation) — The relation p a (T,p) 
as defined in the previous section is multivalued in terms of the density vari- 
able p (fig. 4.1(a)). Such multiplicity, of course, is not physically realis- 
tic. It is therefore convenient to define a new function p = Pp(T,p) that is 
single valued, although discontinuous, and is a physical representation. We 
define this function as 


P p (T,p) 


P a ( T, p) 


= P V ( T ) 


T > T and 0 ^ p ^ 0.0^2 


T < T 


T < T 


and p <; p sv (T) 
or 

p 2 p sl (t) 

and P S y(T) < p ^ P sl (t) 


(4.4) 


where p v (T) is the vapor pressure equation later defined by equation (4.9), 
p sv (T) is the density corresponding .to the saturated vapor boundary, and 
pgL^) is the density corresponding to the saturated liquid boundary. The pro- 
cedure that leads to the evaluation of these quantities is discussed in a later 
section. Note* that although pp(T,p) is identical to p a (T,p) almost every- 
where, it is not convenient to code both these quantities as the same subpro- 
gram. The reason lies in the fact that p a (T,p) is required in the evaluation 
of the implicitly defined boundary quantities pgy and pg L contained in 

Pp(T,p) . 

Caloric equation of state — Once the thermal equation of state and an 
expression for low-density specific heat are given, an equation for the inter- 
nal energy can be derived that is valid for all temperatures above critical and 
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for temperatures below critical provided that the computed pressures do not 
exceed the vapor pressure p v (T) (or the specified density is always less than 
the saturated-vapor density pgy(T)). Introduction of the Clausius -Clapeyron 
equation allows one to generalize the energy equation to include also the low- 
temperature, high-density (i.e., high pressure p > p v (T)) cryogenic region. 

The derivation procedure will not be discussed here (see, e.g., Hust and 
Gosman, ref. 9). The emphasis here is on the presentation of optimal formula- 
tions of the thermodynamic quantities that lead to their efficient computer 
evaluation. 

The caloric equation for internal energy (see, e.g., ref. 9) with terms to 
account for energy changes in the two-phase region is given by 

e = U(T,p) 


= U(T o ,0) + U°( T) - J/°(T 0 ) + t/(T,p) - U( T,0) , T > T 


(4.5a) 


= U(T ot 0) + U°( T) - y°(T n ) + £/(T s p)[S(p sv -p) + S(p-p QT )] 


SL ' 


-y(T,o) + s(p-p sv kz/(T,p sv ) - [/(t,p sl )s(p-p sl ) 


P SL P SV 

P SLi 

( p p sv \ 

P SV P SL 

P 1 

^. P SL -p SV/ 


S(pg L -p) + S( p-pg^) 


dp (t) 

P (T)- T 


dT 


T < T 


c 

(H.5b) 


where the constant U(T o ,0) = Hrp o - RT 0 = 1133.655 J/g-mole (H^ q = 1590.92 9 

J/g-mole and T Q = 55° K) is the low-density or ideal gas reference energy for 
oxygen. The quantity S(x) is the symmetrical Heaviside unit step function, 
which is unity or zero depending on whether the argument x is greater or less 
than zero, respectively , and has value 1/2 when the argument is zero. The 
energy quantity [/(T,p) is obtained by evaluation of the indefinite isothermal 
integral 


u{ T,p) = / P 


P - T 


3T 


PJ 


§£. 


(4.6) 


The ideal gas (or low density) constant-volume specific heat and its indefinite 
integral ^°( T ) (needed in eq. (4.5)) are given by 



(c 9 /T) I 2 


[exp(c 9 /T)-l J 


C 3\ C 8 

+ ( RT ) £n T + R 


(c 9 /T) 


exp(c 9 /T)-lJ 


(4.7) 


(4.8) 
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The coefficients Cg/R are listed for oxygen in the first column of table 4.5* 
Note also that the third term of the second equation has been removed from the 
summation sign. 


Equation (4.5) also requires the vapor pressure equation and its 
derivative as given by 


to p (T) = a + l (a T)* 
v ° tol * 

^ <T) - P * (T> i 


dT 


(4.9) 

(4.10) 


SL=l 


The coefficients are listed in the second column of table 4.5* 


TABLE 4.5.- DIMENSIONLESS COEFFICIENTS FOR THE LOW-DENSITY 
SPECIFIC HEAT RELATIONS (EQS. (4.7) AND (4.8)) AND FOR 
THE VAPOR PRESSURE RELATIONS (EQS. (4.9) AND (4.10)) FOR 
OXYGEN 


i 

c £ /r 

a i 

1 

-1.86442361D 02 

-6.25967185D 01 

2 

2.0784o24xd 01 

2.^7^50^290 00 

3 

-3. 1+261+29 11D-01 

-4.68973315D-02 

4 

3.50297163D 00 

5.48202337D-04 

5 

2.05866482D-07 

-4.09349868D-06 

6 

-1. 110357990-08 

1.91471914D-08 

7 

2.0 8612 876D-11 

-5.13113688D-11 

8 

9 

1.0189U691D 00 

2.23918105D 03 

6.02656934D-14 


(Coefficients from ref. 3) 


The final quantities are the temperature-dependent saturation densities 
p sv (T) and pg L (T) described earlier. 


The caloric equations (4.5) are actually more general than the representa- 
tion given by Hust and Gosman (ref. 9 5 eq. (43)). Here the equation is gener- 
alized to account for the changes in internal energy in the two-phase region 
due to the work done as isothermal changes in volume occur (i.e., the region 
T < T c and pgy < p < Pgp). The factor 



P ) + S(p 


P 


SL 


) 


was introduced to take account of these effects. Because of the coefficient 
step function S(p - PgyK the curley bracket expression makes no contribution 




A 
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to the internal energy for densities less than the value on the vapor "boundary 
p sv (T). For densities larger than that of the saturated vapor the bracket fac- 
tor increases in value as the ratio of specific volume differences given by 


Ay(p >P s V ) 1/Pgy - 1/P _ PsL (p - p Sv| 

Av ^ P SL’ P SV^ 1//p SV _ 1//p SL P ^ P SL“ P sW 


(4.11) 


Equation (4.1l) increases varying linearly with volume until the saturated- 
liquid boundary is reached. The bracket expression has the value unity for all 
values of density in the liquid region p > p SL . The coefficient of the term. 
£/(T,p) in equation (4.5b) is nonzero only outside the two-phase region. Aside 
from the differences in notation, the caloric equation gives the same results 
as that given by Hust and Gosman everywhere outside the two-phase region. 

The Heaviside symmetrical unit step function was introduced to minimize 
the number of regions to be considered by separate energy equations. When the 
energy equations are coded for electronic computer evaluation, however, the 
unit step function is most conveniently represented by use of logical state- 
ments in such a manner that the respective coefficient quantities are evaluated 
only when the step function is not zero (i.e., has positive arguments). 

The quantity in the caloric equation that depends explicitly on the ther- 
mal equation of state is U( T,p) as defined by equation (4.6). It is worth- 
while to obtain an expression for this quantity that depends only on the coef- 
ficients contained within the modified virial equation of state. We first 
write the expression for the thermodynamic derivative (3p/3T) . We find 


where 




+ Z 


8 

l 

i=l 


8 

I 

i=l 


[TA! l) (T)]g! l) (p) 
1 1 

A^ 2) (T)gj l) (p) 


(4.12) 

(4.13a) 

(4.13b) 


( 1 ) » 

The coefficients TA^ (T) contained in equation (4.13a) (prime denotes dif- 
ferentiation) are listed in the second column of table 4.3. Note that these 
coefficients contain terms that differ by integer factors from the terms con- 
tained in the thermal equation of state listed in the first column (except the 
eighth coefficient). Hence, little additional computer time need be expended 
provided the caloric and thermal equations of state are both evaluated together. 
A similar comment applies to the evaluation of the vapor pressure equation and 
its associated derivative, as well as the low-density specific heat and its 
associated indefinite integral. The indefinite integral, equation (4.6), can 
now be evaluated, with the result 
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(4.i4) 


where 


U(T,p) = -RT l Af 2) (T)g! 2) (p) 
i=l 



<p) = f 





The quantities 


gi 2, <p) 


are defined in the second column of table 

( 2 ) ( 2 ) 


that both coefficients Aj_ and ' have been defined such that 
tive products are also dimensionless. This turns out to be parti 
advantageous when units conversion is required. 


4.4. Note 

their respec- 
cular Tv- 


Quantities Defined Implicitly 

Often a choice of independent variables other than density and temperature 
( e . g . , pressure and temperature or energy and density) would be more advanta- 
geous. Since the state equations are sufficiently complicated to preclude the 
existence of explicit relations with the desired variable dependence, one has 
to resort to numerical evaluation. A procedure is described here for "invert- 
ing" the thermal equation of state to obtain, in effect, the state relation 
p = p (p ,T) . Extension of the procedure to obtain analogous equations with 
still different variable dependence is straightforward and is not discussed 
here. 


Density is obtained as a function of pressure, and temperature by solving 
for the zeros of the equation 

f(p;T) = p - P a (T,p) (4.16) 

This equation is readily solved for all values of pressure and temperature, 
except that an additional complication results in the case where temperatures 
below critical are specified. The "successive linear interpolation procedure" 
described in reference 7 was found satisfactory for finding the roots of equa- 
tion (4.l6), particularly since it requires only that the lower and upper 
bounds of the roots be specified. The procedure converges quite rapidly; in 
general, fewer than eight separate evaluations of f(p), or "iterations," are 
required to obtain a value of the root. For specified temperatures greater 
than critical the root of equation (4.l6) lies in the interval 0 < p < 0.042 
g-mole-cm 3 . The constants that represent the interval boundaries were found by 
searching the tables given by Stewart (ref. 3) to find the maximum and minimum 
values of the density. 

In the case that temperatures below critical are specified, f(p) has zero 
value for three distinct values of density. The root-finding procedure then 
requires isolation of the three branches of the multivalued function p a (T,p) 
for T < T c . Figure 4.1 (a) shows the multivalued character of pressure plotted 
as a function of density. The various curves pictured are isotherms. The 
entire region above the critical isotherm is the region where the density 

74 


i 



■bounds have already been specified* Note that a specie is characterized as 
being in the gas phase when stored at pressures and densities that lie on iso- 
therms above the critical isotherm. The region below the critical isotherm and 
to the left of the line labeled P^(p) in figure H.l(a) isolates what we. shall 
call the vapor region or branch. In this region, the pressure is also less 
than the vapor pressure py(T). We define Pg(p) hy 

p B (p) = p l (i^) a (4 ‘ 1T) 

where 


a = S,n(p c /p 1 )/to(p c /p 1 ) 
Pi = 1 atm 

Pl = 3*1CT 4 g-mole/cm 3 


This relation is readily found by observing that it represents a straight line 
on a log-log graph between the pair of points having coordinates (pi,pi) and 
(p c ,p c ). The line pg(p) , although somewhat arbitrarily defined, is a simple 
representation that can be rapidly evaluated numerically. 


The two-phase branch is bounded on the left by the line p-g(p) and on the 
right by the three vertical segments p = p c , 150° < T < T c ; p = 0.02, 


120° < T < 150°; and p = 0.03, T < 120° K. In this region, pressures computed 
from the equation p a (T,p) do not agree with experimental data; hence, this 
region is ignored in the iteration procedure. In fact, the pressure depends 
only on temperature in the two-phase region and therefore is independent of 
density. Finally, the liquid branch is isolated by observing that it is 
bounded by the three-segment curve, the critical isotherm and the already spec- 
ified largest density value P^x ~ 0.01+2 g-mole/cm 3 . In this region , the 
pressure is characterized as having values that exceed the vapor pressure 

V T) - 

The root-finding procedure converges most rapidly for bounds that most 
closely surround the value of the root. The values given above are somewhat 
arbitrary, and one can, perhaps, introduce even better bounds. Those specified 
above, however, were found satisfactory. If, as recommended earlier, the 

A-j^^T) quantities (eq. (1+.2)) are reevaluated only when the temperatures are 

different (a similar comment applies for the gj^Cp) relations), then, 
although about six to eight iterations are required to find p(p,T), the actual 
computer time expended is only a factor of two or three greater than the long- 
est evaluation time required to find p a (T,p) (i.e., the time required when 
both p and T have values that are different when compared to a previous 
evaluation) . 
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For the reader’s convenience the diagram below shows a logical sequence of 
operations that leads directly to the categorization of a specified combination 
of the variables p and T . into a particular branch, . The rounded blocks 
labeled f(p) (eq. (4.l6)) denote the root-finding procedure described in refer- 
ence 7 with appropriate bounds as specified. 

By slightly altering the proce- 
dure just described one can also 
obtain the saturation densities 
p S y(T) and Pgx/T). As commented ear- 
lier, these quantities are the exact 
boundaries of the two-phase region. 
Specifically, these quantities are the 
outer two of the three roots that sat- 
isfy f(p) = 0 , T < T c where the 
specified pressure is that computed 
from the vapor-pressure equation 
p = py ( T ) . We have been concerned 
with the removal of the multivalued 
character of p a (T,p) so that p(p,T) 
can be evaluated unambiguously. It 
turns out that special allowance must 
be made for the computation of the 
saturation densities. This is indi- 
cated on the flow chart by two special 
entry locations (the small circles 
Flow chart that shows the sequence of operations denoted separately as pgy(T) and 

required to evaluate density p(pT) as well as the Psi/-^ ) at ^bich point p = py(T) and 

saturation densities p S y(T) and pg L (T). T are specified. By this procedure, 

the saturation densities can be accu- 
rately computed for all temperatures 

less than about 150° K. For temperatures larger than this value but less than 
T c the saturation-density curves (fig. U.l(b)) vary little with changes in 
density; therefore accurate numerical evaluation is difficult. For the temper- 
atures 150 ^ T ^ T q , these quantities are specified explicitly, a procedure 
that was also found necessary by Weber (ref. 4). The present procedure devi- 
ates from that of Weber in that to maintain internal consistency, special 
tables, with relatively few numbers of elements, were generated by smoothing 
the pgy(T) and Psi/ T ) appropriately obtained from solutions of equation (^.l6) 
(the uppermost value of density in both tables is precisely the critical value 
p c ). A table-lookup procedure with linear interpolation was then used to find 
such values of the saturation densities. 


Quantities Based on Spatially Dependent Variables 

When a system is not in thermodynamic equilibrium (see, e.g., ref. 10, p. 
60 ) much insight on its properties can be gained by the calculation of certain 
thermodynamic averages, such as density and energy. The final state can often 
be predicted if some property of the system (say, energy) is conserved while 
relaxation to equilibrium occurs. In this section we introduce a method 
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whereby one can rapidly predict the final state of a stratified fluid should 
the fluid be subjected to some mechanism that results in complete and adiabatic 
mixing* Since the procedure requires that evaluations be done numerically, 
variables (and summations) defined at discrete points rather than continuous 
variables (and integrations) will be considered. Also, since considerable 
amounts of computer time can be expended in the evaluation of the thermodynamic 
relations over a field of many points , particular emphasis is laid on the 
development of a procedure that minimizes redundant arithmetic operations and 
is therefore relatively efficient. 

The problem described above can be briefly outlined as follows. One first 
finds the quantities 


p(x,y) = p[p, T(x ,y ) ] (4.18) 

e(x,y) = U[t(x , y ) , p(x,y)] (4.19) 

by the use of procedures described earlier. The notation shows results on the 
left-hand side of the equation that are found from the function evaluations 
designated on the right-hand side. The pressure p and the temperature 
T(x,y) are considered known. The total mass M and the total energy E are 
found next by evaluating the integrals 

M = / p dv = pv (4.20) 

T 

E = / pe dv = pev (4.21) 

T 

over the entire tank volume. Introduction of the volume averages p and pe 
for the mass and energy densities as defined above (v^ is the total tank 
volume) then allows one to find 


T cq1 = T(pe/p , p) 

p col = P p (T col’ ?) 

The temperature is found by inverting the caloric equation of state (4.19). 

The results denoted by the symbols T co j_ and p co l are the temperature and 
pressure to which the system would "collapse" and that uniquely characterize 
the uniform final state of the oxygen contained within the tank. During a 
space flight, p remains constant unless fluid is drawn from the tank. If 
severe thermal stratification occurs, mixing (by vehicle maneuvers, for exam- 
ple) can bring about large abrupt pressure decreases (see fig. 4.1(c) and dis- 
cussion) that could conceivably lead to spacecraft equipment failures. These 
concepts are not new, however, and are discussed elsewhere (ref. 2 and ch. 6). 


(4.22) 

(4.23) 
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The evaluation of equations (4.18) and (4.19) is particularly time consum- 
ing . The matrix 


where 


Tj k = T(j Ax, k Ay) 


x = j Ax , 1 < j < JM 

y = k Ay , 1 < k < KM 


(4.24) 


(4.25) 


can contain many elements. Furthermore, it can occur that the temperature 
field describes situations where many volume elements are at temperatures dif- 
ferent from that of the uniform surrounding gas. Once the density and energy, 
equations (4.18) and (4.19), are evaluated, we could efficiently evaluate the 
integrals, equations (4.20) and (4.21), if like temperatures were grouped 
together. However, we must consider that the are randomly mixed and like 

temperatures are not ordered in any systematic manner. It also appears waste- 
ful to recalculate values for, say, density and energy when differences result 
that are not numerically significant (e.g. , less than 0.1 percent). In the 
interest of minimizing the number of separate evaluations of thermodynamic 
functions a procedure is to be described whereby, in essence, one first counts 
the number of times that a particular temperature appears in the array T^. 

In this manner a "temperature distribution function" is then constructed that 
applies strictly to an equally spaced Cartesian coordinate system. In non- 
Cartesian systems, however, additional complexity results. The analytical 
foundation for the temperature distribution function is established in the 
discussion that follows. 

We first define a dimensionless temperature 


T = A (T - T . ) 
o mm 

or, in the case that the spatial dependence is discrete. 


T„ = A (T - T . ) 

jk o jk mm 


(4.25a) 


(4.25b) 


where A Q and '•'-min are both constants ; the value of A Q is arbitrary and 
T m in is the very lowest bound of all possible temperature values T (or T^)* 
A temperature distribution function F(t) is now defined such that F(x)dx 
is proportional to the volume of a gas with dimensionless temperature between 
t and x + dx. Then 


F(x) = C / 6(x-x)dv (4.26) 

V T 

where C is an arbitrary constant set equal to unity; the constant can be used 
to normalize the distribution function in such a manner that its Integral over 
all possible temperatures is unity. The function 5(x - x) is the Dirac 6 
function defined such that 


k 
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/ S(t - x)dx 
T 1 


equals unity provided ti < t < x 2 and equals zero if x is not within the 
integration interval. 

If the volume elements are discretized, we have 


F ( x ) dx 


JM,KM 

I [S(t 

j ,k=l 


T jk> - S(t ' T jk + dT)]M Jk 


JM,KM 

l «<T 

j »k=l 


T Jk )dTW Jlc 


(4.27a) 

(4.27b) 


where S denotes the symmetrical Heaviside step function described earlier; 

- Av^N^/V^ a "weighting" function; Avj^. is the gas volume associated 

with the discrete volume element of the cell labeled j ,k ; N-^ corresponds to 
the total number of volume centered computational grid points ; and Vrp already 


defined, is the tank volume. It follows that when C = 1, ^^Wjk “ ^t * To 

further clarify the meaning of the weighting function W^ k) if all Av.^. were 
equal (such that we have an equally spaced Cartesian system), the would 

then all be unity; that is, V T /AVj k = 


We now seek a discretized version of F(t). 
define a function 


A simple possibility is to 


N + f 


F n = / x F(x)dx 


N- f 


JM,KM 

i 

j »k=l 


S|t Jk 





(It. 28) 


Then F^ is the weighted number of computational grid points with dimension- 
less temperature between N - 1/2 and N + 1/2. The mean dimensionless tempera- 
ture x within this interval is N (an integer). How one can evaluate F^ 
rapidly with a computer is explained later. Other thermodynamic variables such 
as the associated dimensional temperature Tjj, density p^, and energy e^ can 
be evaluated according to the relations 


T = T + 


N 


P N 


p(p»V 


S N = 


U( W 


(4.29a) 

(4.29b) 

(4.29c) 
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It will be shown that fewer evaluations of these quantities over the entire 
grid occur than if a procedure is followed whereby the are used directly. 

Consistent with our definition for the distribution function, equation 
( 4 . 26 ), it can be shown that, in general, the volume average of a thermodynamic 
quantity, say g(*r) 9 is given by 


= iT Jq F(T)g( T )dT 

t 


(4.30) 


We separate the integration interval into a sum of subintervals that cover the 
range of possible temperature values 


g = 


w 

1N max 1 

1 v i*N+ - 


N. 


If 2 g( t )f( r ) dr 

_ 1 1\T 1 


T N=1 N- f 


where 


(It. 31a) 


K = A (T - T . ) 
max o max mm 


(4.31b) 


To zeroth order, equation (4.31a) is approximated by 


N 

j_ max 

® ~ nT E 

t N=1 ^ 


(4.32) 


The accuracy of the average computed by this equation can be improved by 
increasing the value of the arbitrary constant A Q , but at the expense of 
increasing the number of temperatures % at which the therodynamic functions 
are computed. 

A more efficient procedure is to use a higher order method as follows. 

The function g(x) can be expanded in a Taylor series 

g(r) = l & ( T ~ (4.33) 

u m! 

m=o 


Substitution of the series into the exact equation (4.31a) yields 


where 



N 

max 00 


1 1 


s (m) W f 

m! Nm 



- N) m dr 


(4.34a) 


(4.34b) 
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Substitution of equation (4.27b) into equation (4.34b) leads to 



I 






(4.35) 


With equations (4.34a) and (4.35) > it is possible, in principle, to compute g 
to any order of accuracy, although if one carries too many terms, the original 
purpose for introducing the procedure - increased computational efficiency - is 
defeated. With m = 0 , we see that equations (4.34) and (4.35) become identi- 
cal to equation (4.28). One can readily evaluate an equivalent set of rela- 
tions that are valid to first order. One first introduces a simple difference 
relation for the derivative g( 1 )(N) and then with the definition 


F = / N+1 F(t ) ( l - | t-N | ) dx (4.36a) 

N-l 

it can then be shown that the average 

N 

max 

g = I F g (4.36b) 

N=1 

is correct to order (t-N) in the expansion (4.33) for g(x). 

The temperature distribution functions are particularly simple to find 
with an electronic computer especially when integer subscript notation is used. 
To illustrate this point, let F(N) be the variable in FORTRAN language that 
represents the zeroth-order representation for the distribution function, equa- 
tion (4.28). The one-dimensional array F(N) must be of sufficiently large 
dimension to contain all of its possible elements; it must contain at least 
NMAX elements, where NMAX is given by equation (4.31b). We next set all of 
the elements of F(N) equal to zero since those locations in computer memory 
will be used to accumulate the sums denoted on the right-hand side of equation 
(4.28). We then process each element in the array T(j,K) that represents the 
variable T jk* We d -° this by first evaluating N = A 0 [T(j,K) - TMIN] (note 
that we are using the automatic truncation feature of "FORTRAN language whereby 
the integer on the left-hand side of the equation represents the largest inte- 
ger value less than or equal to the floating-point result of the expression 
given on the right-hand side of the equation). We define TMIN as the FORTRAN 
equivalent of T m j_ n . The weighting factor W(j,K) associated with the spatial 
point j ,k is then added to what is contained in F(N) by the use of the FOR- 
TRAN statement F(N) = F(N) + W(j,K). This procedure, briefly outlined, is 
easily generalized such that one can find the associated distribution function, 
equation (4.36) or the higher order forms, equation (4.35)* Several example 
plots of temperature distributions are given in chapter 6. The procedure for 
finding the averages, equation (4.32), or, in the case of the higher order rep- 
resentations, equation (4.34a), is straightforward and simple to program on an 
electronic computer. It need only be said that one will also obtain additional 
savings in computer time if the thermodynamic quantities, say p^, e^ , that 
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correspond to the temperatures Tjj are not computed when the associated dis- 
tribution function quantities F(N) have zero value. 


EXAMPLE EVALUATIONS OF THE THERMODYNAMIC FUNCTIONS 


The explicit representations p = p a (T,p), p = Pp(T,p)., and e = U(T,p.) 
have been defined. By the procedure discussed, one can invert the above rela- 
tions so that, in effect, the above list is readily expandable to also include 
representations of the type p = p(p,T), p = p(e,T), T = T(p,p), and T = T(e,p). 
Our number of equations is now sufficiently complete that a value of any one 
variable of the set p, p, T, and e can be found once values are specified for 
any other two variables of the set. In principle, by introducing additional 
thermodynamic functions (such as enthalpy and entropy, for example) the list 
can be expanded further. 

In this section a series of plots are given that illustrate the flexibil- 
ity of the procedures discussed. These plots are arranged so that the curves 
within each graph represent cuts along curves on both the p-p-T (figs. 4.1, 

4.2, and 4.3) and the internal energy surfaces (figs. 4.4, 4.5, and 4.6). The 
figures are further arranged according to the particular variables displayed 
along the ordinate and the abscissa. Pressure is plotted versus density in 
figure 4.1. Other figures show pressure versus temperature (fig, 4.2), density 
versus temperature (fig. 4.3), and energy versus density (fig. 4.4), tempera- 
ture (fig. 4.5), and pressure (fig. 4.6). Each figure shows curves for con- 
stant temperature (isotherms), constant pressure (isobars), constant density 
(isochores), or constant energy (isoenergetics), as appropriate. Such plots 
may be instructive as well as useful as a ready reference. 

Recall, however, that the basic equations used are semi empirical and 
therefore the results plotted have quantitative value only within the regions 
where the constants were fitted by the least-squares procedure with experimen- 
tal data by Stewart. Thus, the extreme values that are displayed (in particu- 
lar, for the largest values of temperature and density) may not be precisely 
accurate. The less accurate ranges of the curves are often easily identified 
by the inflections or by the discontinuities in curvature occurring at the very 
ends of the curves. In the case that more accurate values are required than 
can be read from the curves, the reader is referred to the tables given in 
reference 3 or 4, or the equations given in this chapter. 

Figure 4.1(a) shows a number of isotherms on a log-log plot of pressure 
versus density that illustrate the analytical behavior of Stewart's equation of 
state p a (T,p) ((4.1) , or, as modified, eq. (4.2)). The lowest isotherm corre- 
sponds to a temperature of 100° K. The higher curves are for higher tempera- 
tures that increase by increments of 10° K and by 20° K for temperatures above 
200° K, except the curve drawn for the critical isotherm T = T c , between the 
isotherms 150° and l60° K. Small crosses are drawn on the isotherms below the 
critical to show where the phase boundaries intersect. Knowledge of the exact 
behavior of these isotherms as computed from Stewart's equation can be quite 
helpful. Since the smooth multivalued character of the equation is displayed. 


one can readily isolate appropriate "branches in the manner indicated by the 
straight-line segments so that an inverse of the function can he obtained, as 
discussed in the previous section. Although the state equation contains trans- 
cendental terms of a complicated nature fsee, e.g., eq. (4.l)) their plotted 
results have the same overall characteristics of an equation that is a cubic in 
density for the range of parameters shown. As such, the equation is not 
greatly different from the van der Waals equation (ref. 6) except, of course, 
for the magnitude of the computed values. One observes that the critical iso- 
therm contains an inflection at the critical point (which can be identified as 
the point where the line segments converge) as indeed it must. An inflection 
also appears at the critical point in the density-temperature plane shown by 
figure 4.3. 

Figures 4.1(b), 4.2(a), and 4.3(a) illustrate the behavior of the "physi- 
cal” representation of the thermal equation of state (4.4). Note that pressure 
is independent of density in the two-phase region as it should be. The manner 
in which the two-phase region degenerates to the vapor pressure curve can be 
understood by comparing figures 4.1(b) and 4.2(a) (also, cf. fig. 4.3(a) with 
fig. 2 in Weber T s paper,' ref. 4). Difficulty will be encountered if one tests 
the computed pressure p = p a (T,p) against the vapor pressure py(T) as a pos- 
sibly simpler scheme than that described with regard to equation (4.4) for 
separation of the vapor from the liquid regions (the rule quoted in elementary 
texts on thermodynamics, such as ref. 5* is that for T < T c one has vapor for 
p < py ( T ) and liquid for p > py(T)). The test fails in the comparison of 
p^(T,p) with p-y(T) since the pressure computed from the analytical representa- 
tion has widely varying values in the two-phase region. One must use the com- 
puted saturation densities (see flow chart) to correctly separate regions. 

Figures 4.1(c), 4.2(b), and 4.3(b) show the behavior of isoenergetic 
curves on the respective coordinate planes p-p, p-T, and p-T. In comparing 
figures 4.1(b) and 4.1(c) at the lower densities, note that as the temperature 
increases, the effect of density on energy becomes less (the isoenergetic and 
isothermal lines become colinear implying e *> e(T) = constant). We recall 
from elementary principles (see, e.g., ref. 5 or 6) that in the low-density 
gaseous region the energy depends only on temperature e = e(T). This feature 
can also be observed in some of the other figures; see, for example, figures 
4.3(b) , 4.4(a) , and, in particular, 4.5(a) where we see that energy remains 
approximately constant for the isochores less than 1 g-mole/cm 3 . Thus an 
appropriate limiting behavior is correctly displayed. 

One other feature of importance displayed in figure 4.1(c) is that in the 
cryogenic region where pressure and density both have large values, adiabatic 
mixing of a thermally stratified fluid contained with a closed container (here 
one follows a constant energy curve, say, the curve labeled -2000 J/g-mole) 
would have only a small effect on the value of density, but could lead to 
rather severe pressure reductions. 2 


2 By use of thermodynamic arguments it can be shown that during mixing the 

pressure reduces rather than increases (ref. 2). 
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Figures 4.4, 4.5 9 and 4.6 illustrate the effect of isotherms, isobars, and 
isochores drawn on an energy surface. Note that discontinuities in energy 
occur along isotherms in the energy-pressure plane (or along isobars in the 
energy-temperature plane); for example, see figures 4.5(b) and 4.6(b), The 
magnitude of the energy jump is a measure of the heat of vaporization associ- 
ated with the respective values of temperature and pressure noted. Analogous 
discontinuities are not observed in the case of energy versus density 
(fig. 4.4) since the effect of work done in compressing the vapor to obtain a 
liquid was accounted for in the caloric equation. One can use the abrupt 
changes observed in slope values to mark the phase-change boundaries. In the 
case of figure 4.5(a), the slope of the isochores (3U/3T)p is identified as 
the constant-volume specific heat. Plots of specific heat are given in 
reference 3 (pp. 47 and 48). 
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(a) Isotherms computed on the basis of the analytical representation p a (T,p), equation (4.1) or (4.2). 

Figure 4.1.— Pressure versus density. 
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(c) Isoenergetics computed from p(e ,p) obtained by inverting equations (4.5). 
Figure 4.1 Concluded. 
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(a) Isochores. 

Figure 4.2 — Pressure versus temperature. 
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(a) Isobars. 

Figure 4.3 — Density versus temperature. 
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(b) Isoenergetics. 
Figure 4.3.— Concluded. 
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(a) Isochores. 

Figure 4.5 .— Internal energy versus temperature. 
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Figure 4.5.- Concluded. 
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Figure 4.6 — Internal energy versus pressure. 


96 




Energy, J/mole 



Pressure, atm 


(b) Isotherms. 
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5. ANALYSIS OF THERMODYNAMIC STATES RESULTING FROM 


SMALL-DENSITY-VARIATION APPROXIMATION OF 
FLUID MOTION IN CRYOGENIC OXYGEN TANKS 


Barrett S. Baldwin, Walter A. Reinhardt, and Yvonne S. Sheaf fer 



SUMMARY 


Methods are developed for analysis of the thermodynamic states predicted 
by a small-density-variation approximation of the fluid motion in a cryogenic 
storage tank. Relations are derived for evaluating variations in pressure and 
other thermodynamic quantities resulting from conduction and convection of 
heat. Methods are developed for simulating the heater in an Apollo oxygen tank 
in a manner compatible with the small-density-variation approximation of the 
fluid motion. Correction procedures are specified to account for compressibil- 
ity effects and to suppress nonphysical behavior introduced by the numerical 
method used to determine the conduction and convection of heat. 


INTRODUCTION 


The integration procedure for evaluating the motion and temperature dis- 
tribution of the fluid described in chapters 2 and 3 is based on several 
approximations. Although pressure gradients must be considered in the momentum 
equations, it can be shown that for the low- velocity cryogenic oxygen flows 
under consideration, pressure gradients can be neglected in the energy equa- 
tion. Similarly, the variations in pressure with time do not affect the motion 
of the gas directly, but have only a cumulative effect in the energy equation. 
The most important coupling between the momentum and energy equations is 
through temperature gradients and convection of temperature variations. As a 
result , for the purpose of finding the motion of the gas , the density can be 
considered a function primarily of temperature, with pressure as a slowly vary- 
ing parameter. If the temperature variations are sufficiently mild, the den- 
sity and enthalpy can be approximated by the first several terms of series 
expansions in temperature with coefficients that vary slowly in time due to 
cumulative changes in pressure and temperature level. 

In chapter 2 it is shown that to lowest order the flow equations can be 
cast in a form that is independent of pressure and density variations. The 
associated, integration procedure described in chapter 3 requires as input a 
mean density p and thermal expansion coefficient 6, each of which can be 
time-dependent to allow for cumulative changes in pressure and density level. 

A time- dependent temperature distribution results from this integration pro- 
cedure. In this chapter, we indicate how this output can be used to determine 
the attendant slowly varying pressure as well as other thermodynamic variables 
such as the density distribution and the potential pressure decay that would 
result from adiabatic mixing. Methods for simulating the heater in an Apollo 
oxygen tank are developed. Correction procedures are also developed to account 
for compressibility effects and to counteract nonphysical behavior introduced 
by the numerical method (ch. 3) used to compute evolutions of temperature. 

So far as the thermodynamic state of the gas is concerned, pressure gradi- 
ents are negligible at the low velocities that occur in a cryogenic oxygen 
storage tank in a near zero gravitational environment. The pressure can there- 
for be considered uniform in the tank. The kinetic energy associated with the 
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motion of the gas and the dissipation term in equation (2.97) are also negligi- 
ble compared to changes in the internal energy of interest. Thus the problem 
here is to evaluate the time- dependent thermodynamic state of a stationary 
stratified gas from a knowledge of the time -dependent mean density and tempera- 
ture distribution. The degree of rigor brought to bear is independent of the 
means used to arrive at the input information on mean density and temperature 
distribution. The methods developed for this purpose therefore would be applic- 
able to the results from more accurate higher order or three-dimensional anal- 
yses of the convection and conduction processes within the tank. Regarding the 
rigor of the analysis, we have considered the problem at two levels. In pre- 
liminary work, the van der Waals equations of state were used as described in 
the next section. In subsequent sections, we include results based on the more 
exact thermodynamic relations developed in chapter k . 

APPROXIMATE PROCEDURE BASED ON VAN DER WAALS EQUATIONS 


The problem at hand is to find the pressure and other thermodynamic quan- 
tities when the mean density and temperature variation in the tank are known. 

In this section, we describe the procedure used to obtain preliminary values of 
these thermodynamic properties by an approximate method based on the van der 
Waals equations of state. The mean density p is known in terms of the total 
mass of oxygen M from the relation 


p = M/v t (5.1) 

where Vy is the tank volume. The rate of change of mean density is found by 
differentiation 


dg_ _ dM/dt 

at v T 


where dM/dt is the rate of gas removal (typically 0 to 3 lhm/hr) . Values of 
p are computed at each time step from the relation 


^ P L+At 


, At ( dM/dt ) 
(p) + + w 


(5.3) 


using input values of dM/dt and V^. 

At each time step the values of temperature T^ at the computational 
grid points (see fig. 3.1) are determined from the integration procedure. The 
corresponding densities p ^ are computed to lowest order from the relation 


P = p[l - B(T - T)] 


(5.4) 


where T is the volume average of the temperature variation Tj^ and 
6 = (-l/p)(3p/3T) is the coefficient of thermal expansion corresponding to 
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the state p, T (see eq. (2.1 6 )). The van der Waals equation for pressure is 

(5.5) 


RTp o 

P = H 1 r - ap z 


1 - hp 

Differentiation leads to the relation 

R(1 - hp) 


3 = 


RT - 2ap(l - hp) 


(5.6) 


for the coefficient of thermal expansion. In our preliminary results reported 
at meetings of the Apollo Cryogenic Oxygen Tank Analysis Team at the MSA 
Manned Spacecraft Center, a constant value of g "was used that was evaluated 
at the initial values of p and T. 


Since a uniform spacing of computational grid points is employed in the 
integration procedure , the volume average T can he computed according to the 
relation 


T = 


l 

jW 


I 

j >k 



(5.7) 


where is a weighting function (proportional to the volume associated with 

each grid point) that is taken equal to 1.0 at interior points, 0.5 at boundary 
points except in the corners, and 0.0 in the corners. The corners are excluded 
(zero weight) because in the integration procedure described in chapter 3 tem- 
peratures at the corners are not used and are not computed. The error intro- 
duced by this omission is negligible for the 17 x 17 array of grid points used 
in the calculations . 

Substitution of the Tj^ and into equation (5*5) would lead to pres- 

sures pjk that are not all equal as they should be. Since the pressure 
should be uniform in the tank, it is expedient to evaluate an average pressure 
according to the relation 


(p + ap 2 )(l - hp) = RTp 


where the bar indicates a volume average of the same type as in equation (5*7). 
Substitution of equation (5- 4) into th is relation, use of the fact that 
(T - T) =0, and omission of (T - T) a terms found to he negligible in the cases 
considered lead to the approximate expression 


- = RTp_ _ a( - )2 + [a(p) 2 (3bp - l)g 2 - Rpg] (T - T) 2 
1 - hp 1 - hp 
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for the pressure in terms of average density, average temperature, and mean 
squared deviation of the temperature from the average value. The latter quan- 
tity is computed at each time step as in equation (5.7) according to the 
relation 


(T - T) 2 = 


I 

.1 



- T) 2 W._ 
Ok 


l 

0,k 



(5-9) 


In addition to the pressure, the potential pressure decay (pressure decay 
that would result from complete adiabatic mixing) is of interest. This quan- 
tity provides a conservative estimate of the maximum drop in pressure that 
could result from an abrupt spacecraft maneuver. This pressure drop could 
result in the tank fluid becoming a two-phase mixture, an undesirable state 
according to design criteria. The potential pressure decay that arises from 
nonuniform heating under various operating conditions is thus an important 
indication of the level of stratification of the fluid. In general, the pro- 
cedure for evaluation of this quantity is as follows : 


total mass 

M 

= pVij 

total internal energy 

E 

= epVrji 

specific internal energy of collapsed state 
resulting from complete adiabatic mixing 

e col 

= E/M 

density of collapsed state 

P col 

= P 

temperature of collapsed state 

T col 

= T * e col’ p col 

collapse pressure 

■^col 

= p(T ool’ p ccl 

potential pressure decay 

(p - 

Pcol 1 


where e is the specific internal energy at individual grid points and ep is 
a volume average. The functions T ( e col > p ^) P ^ T col* P col^ are to " be 

found by inversion of the equations of state for specific internal energy 
U(T,p) and pressure p(T,p). 

The above procedure for determination of the potential pressure decay can 
be carried out using approximate equations of state or the more exact relations 
given in chapter k. Our preliminary results were based on the van der Waals 
equations of state, equation (5*5) and (ref. l) 

e = ci + c^T - ap (5-10) 


Substitution of equation (5-5) into (5.10) to eliminate T yields 


c 

e = ci + — (p + ap 2 )(l - bp) - ap 

This relation can be used to evaluate ep* and obtain an expression for e co q 

in terms of p, p, (p 2 ), and (p 3 ) according to the foregoing general procedure 
for finding the potential pressure decay. Direct substitution in the above_ 
relation also yields an expression for e co q in terms of P C ol an( 3- Pool = p * 

Equating the two expressions for e co q then leads to the formula 

ab[^- (p)3] - a[l - (R/c v )][72 - (p) 2 ] 

P ~ P col ~~ 1 - bp 


where p 2 and p 3 indicate volume averages. Substitution of the first-orde r 
density relation equation (5*4), rearrangement, and omission of (T - T) 3 
terms leads to the approximation 


a(p ) 2 [ 3bp - 1 + (R/c v )]3 2 (T - T) 2 
P ~~ P col " 1 - bp 


(5.11) 


which expresses the potential pressure decay in terms of the average density 
and mean squared deviation of temperature from the average value. 

The specific heat at constant pressure is needed in the integration pro- 
cedure described in chapters 2 and 3. The van der Waals equations of state 
(eqs. (5.5) and (5*10)) can be used to derive the relation 

C p = C V + 1 - 2ap(l - bp)2/RT ^• 1 ‘ 

In our preliminary results a constant value of c p is used that is evaluated 
at the initial values of p and T. The values of the constants R, b, a, and 
c v used in the van der Waals equations were chosen such that the critical 
pressure, temperature, and density are matched exactly according to relations 
given by Hirschfelder , Curtiss, and Bird (ref. l). 


PROCEDURE BASED ON EXACT THERMODYNAMICS 


In this section a more exact method is developed for finding the pressure 
and other thermodynamic quantities when the mean density and temperature dis- 
tribution in the tank are known. This method makes use of a temperature dis- 
tribution function to save considerable computation time. A method for simu- 
lating the heater in an Apollo oxygen tank also is described, and correction 
procedures are developed to account for compressibility effects and to counter- 
act nonphysical excursions of pressure introduced by the numerical method 
(ch. 3) used to compute evolutions of the temperature distribution. 
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Computation of Tank Pressure 

When the mean density and temperature distribution are known, the uniform 
pressure in the tank can be computed. For this purpose it was found expedient 
to determine a temperature distribution function F^ of the type developed in 
chapter k (after eq. (U.l8)), rather than considering thermodynamic quantities 
at each computational grid point. The procedure using a distribution function 
and computation of thermodynamic quantities at temperatures Tjj replaces com- 
putation at grid temperatures Tj^. The advantage is that there are many more 
values of Tj k than of Tj$, so many unnecessary repetitions of nearly identical 
calculations are eliminated. The temperature distribution function F is a 
function of temperature to be evaluated at a series of fixed, equally spaced 
temperatures Tjj with H = 1, 2, 3, . . . . N max . ®h e temperature dependence 
of F is here indicated parametrically by expressing F as a function of IT, 
that is, Fjj. The function F$ is defined to be a weighted number of computa- 
tional grid points with temperatures between % - AT and Tjj + AT where the 
Tjy constitute a fixed array of temperatures with uniform spacing equal to AT. 
The weighting employed is proportional to the volume associated with each com- 
putational grid point. Interior points are given a weight W^ k = 1.0, for 
boundary points = 0.5, and the corner points are given zero weight since 

their temperatures are not computed in the integration procedure. For each 
value of j and k, Wj^ is assigned to the two Fjj between which its tempera- 
ture lies in proportion to its proximity to each. That is, if 
% < ^jk < %+b % increased by an amount Wj^Ct^+i - Tj k )/AT, ana % +1 
is increased by an amount Wj fc (Tj k - T^)/AT. Thus the sum 

N 

max 


is equal to the total number of interior computational grid points plus half 
the number of boundary points, not counting corners. The quantity Fj^ is 
essentially equal to the total number of computational grid points with temper- 
ature between % - ( 1/2) AT and Tjy + ( 1/2) AT except for a small readjustment 
corresponding to a linear interpolation. Figure 5-1 shows an example of a tem- 
perature distribution Fjy plotted versus T^. Additional examples and further 
discussion of the meaning of this distribution function are given in chapter 6. 

An array of temperatures T^ is associated with the temperature distribu- 
tion function Fjj. Since the pressure is "uniform in the tank, at a given tank 
pressure p, associated arrays of density and internal energy e^ can be 

computed according to the relations 

P w = (p,T n ) (5.13) 

e N = U(T N’ P N } (5 ' l4) 

The computation of the function p(p,Tpj) is described after equation (4.l6) and 
U(Tjj,pjj) is given by (4.5). Since the % are proportional to the volume of 
gas in the temperature range Tjj - ( 1/2) AT to Tjy + ( 1/2) AT the volume averages 
p and pe can be computed according to 
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max 


pe = 


I 

N=1 

Nmax 

l 

N=1 


N 


F N P N e N 


N 


(5.16) 
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Figure 5.1.— Typical temperature distribution function. 


At the beginning of an integration 
for the motion and temperature of the 
fluid a tank pressure p and values of 
temperature at the computational 

grid points are specified. The distri- 
bution function F^ is computed, and 
the above relations are used to deter- 
mine the initial values of p and pe. 

In the subsequent integration the vari- 
ation of p is computed from equation 
(5.3) and depends on the specified 

rate of gas removal cLM/dt. The values of T- k resulting from the integration 
for the motion of the gas (see eqs . (3. 2d) and (2.95a)) are used to compute F^ 
at the end of each time step. The problem then arises of computing the change 
in pressure Ap in each time step. It should be recalled that pressure p is 
uniform in the tank. By differentiation it can be seen that changes in the p^ 
are related to changes in pressure according to 


. Ap 

Ap N (3p/3p) r 


N 


(5.17) 


since the values of Tjj are held fixed throughout the integration. Substitu- 
tion of this in the relation 


(pN ) n+1 = (PN)„ + a PN 


n 


yields 


(pn) 


(pn)„ + A P 


n+1 (3p/3p), 


l N 


(5.18) 

(5.19) 


where the subscript n refers to the time step in the integration. Multipli- 
cation of the last equation by (Fjj) n+1 , summation over IT, and substitution of 


\ax _ N max 

^^n+l^ PN ^n+i = p n+l ^ ^ F N^ n +i 


( 5 . 20 ) 
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yield 


^max \ax Vax 

VjJ, ( %>n+. W,.) 1 '! 1 .*®,!, (%) n+1 0p/3p) % (5.21) 


Since P n+1 and - ( F N) n+1 are known at the end of the (n+l)th time step, this 

equation can he used to compute the change in pressure Ap associated with the 
time step. Once Ap is known, equation (5-19) can he used to compute the 
updated values of density a c keck, P n+1 can he computed using 

equation ( 5 . 20) for comparison with the imposed value P n+1 from equation (5*3). 
Finally, the updated pressure is computed according to 

Pn+1 = p „ + 4p (5.22) 


The structure of equation (5.21) is such that the value of P n +i computed 
from equation (5.20) will always he driven toward the imposed value from equa- 
tion (5.3) so that cumulative drifts cannot occur. This follows from the fact 
that the first term on the right of equation (5*2l) closely resembles the left 

side of equation (5.20). Cumulative drifts of the individual (pm*) computed 

^ n+i 

from equation (5*19) can take place, however. To avoid this, at every tenth 
time step, the are recomputed according to equation (5.13). It has been 

found that the foregoing procedure is quite stable, and, for the sizes of time 
step imposed by stability criteria of the integration procedure, it is quite 
accurate with a 1° K spacing of the temperature elements T^. 
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Once the pressure is determined the potential pressure decay can be com- 
puted by the procedure described in chapter k that utilizes the distribution 
function As a check the potential pressure decay can also be computed by 

the alternative method of chapter 4 for which the pressure p and temperatures 

at computational grid points are uti- 
lized. Figure 5.2 shows a comparison 
of potential pressure decays computed 
by the two methods for a linear temper- 
ature distribution. The dashed lines 
indicate the potential pressure decay 
for each case computed exactly by eval- 
uation of the thermodynamic quantities 
at all computational grid points . The 
symbols indicate values of the poten- 
tial pressure decay computed using the 
distribution function F^ for various 
values of the temperature array spacing 
AT. It can be seen from the results in 
figure 5.2 that the temperature array 
spacing AT = 1° K provides adequate 
accuracy for all the cases considered. 
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Figure 


5.2.— Effect of temperature array spacing on 
computed potential pressure decay. 
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Method for Heater Simulation 


An internal heater is used in the Apollo oxygen tanks to increase the 
pressure when it falls below 865 psia due to gas removal. In unpublished work 
of C. K. Forester, D. D. Rule, and H. W. Patterson reported at meetings of the 
Apollo Cryogenic Oxygen Tank Analysis Team, a segment of the wall was used to 
simulate such a heater. It was found that the boundary-layer flow in the 
neighborhood of the heater cannot be adequately resolved with a uniform grid 
spacing when the acceleration field is of order 10" 6 g or greater. Methods can 
nevertheless be found that lead to physically reasonable and qualitatively cor- 
rect results if attention is confined to energy conservation, and accurate 
values of heater temperature are not required. 

In this section, a method for simulating the heater by assigning selected 
computational grid points as heater elements is described. When the heater is 
on, the temperatures of the grid points assigned as heaters are increased by an 
amount dependent on the heat capacity of the gas associated with each grid 
point. The heat capacity depends on the volume of gas associated with a 

grid point, which is equal to the depth of the tank times the product of grid 
point spacings in the x and y directions. Adjustments in heater temperatures 
are made at the beginning and end of a time step in the integration for motion 
and temperature of the fluid described in chapter 3. The energy balance will 
be properly maintained if at the beginning of each time step the heater element 
temperatures are increased by an amount AT h corresponding to a specified 
heater power dQ/dt added to the total gas volume associated with the 

heater element computational grid points; that is. 


before the integration step, 
computational grid point is 


. (as/g>t ( 5 . 23 ) 

E CC p V HT 

The gas volume associated with each interior 


V H " 


V T Ax Ay 

n 

x y 


( 5 . 2 U) 


where Ax and Ay are the distances between grid points ; the depth of the tank 
is equal to Vrp/£ x 2y; Viji is the total tank volume; and are ^h e tank 

dimensions. 

If the heater element computational grid points are interior points, the 
total heater gas volume Vjpp is equal to the number of heater elements times 
Vjj . But heater volume elements on the boundary are half the value given in 
equation (5.2*0 since the boundary passes through the grid points. When the 
heater elements are in the interior their temperatures will rise until a bal- 
ance is reached between the temperature increase from equation (5-23) and the 
decrease due to convection and conduction computed in the integration procedure 
of chapter 3. When the heater elements are on the boundary, however, the inte- 
gration procedure does not modify their values. In that case, it is necessary 
to allow for a decrease in temperature at the end of each time step according 
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to the amount of heat transferred to the interior grid points from the heater 
elements. For a heater on a vail parallel to the y axis the appropriate 
change in vail heater element temperatures is 


4T h 


2k(T„ 


- T 

interior 
pc" (Ax)2 
P 


) At 


(5.25) 


after the integration step. When the heater is turned off (dQ/dt = 0) the 
increase ATg computed in equation (5.23) is zero. If the heater is on a 
vail, the ATjj computed in equation (5.25) corresponds to an insulated vail 
boundary condition vhen the heater is off except for a small (physically cor- 
rect) lag due to the heat capacity of the gas adjacent to the vail. 


The same type of computation as that in equations (5*23) and (5.25) can be 
used at all boundary points to simulate the heat leak from the exterior of the 
tank. In that case dQ/dt in equation (5.23) is replaced by the heat leak 
rate dQ^/dt. 


The pover radiated from the heater can be alloved for by means of the 
relation 


_ o 

dt dt ^rad 


vhere dQ/dt is the specified input heater pover and 


(5.26) 


Vd - e V t h - <5.27) 

vhere e is the emissivity (typically 0.32), a r the Stefan-Boltzmann constant, 
and Ajj the heater area. According to reference 2, less than 10 percent of 
the radiated pover is absorbed in the oxygen (usually much less depending on 
the heater temperature). Therefore, most of the radiation is absorbed in the 
tank vail, and Q ra & should be added to the heat rate dQ^/dt. That is, 
radiation from the heater is absorbed in the vail and heats the fluid adjacent 
to the vail, by conduction in the same manner as heat leaking to the interior of 
the insulated vail from the warmer outer vail of the tank. 


CORRECTION FOR SPURIOUS INTERNAL SOURCES 


The integration method for the motion and temperature of the fluid 
described in chapter 3 utilizes the so-called "conservative form" of the con- 
vective terms. The method is not strictly conservative in the sense of ref- 
erence 3, although it may provide more accurate treatment of convection effects. 
Also the treatment of heater and vail computational grid points described in 
the preceding section is not strictly conservative. For that reason small 
spurious internal sources may lead to erroneous cumulative computed pressure 
variations. To avoid this the method can be made conservative overall by 
applying a correction to the temperatures Tj^ at the end of each time step. 
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For that purpose integration of the energy equation over the tank volume, use 
of Gauss f theorem, and rearrangement lead to the expression 


where 


dE dM 

dt dt o dt 


E = peV T 
h Q s h(p,T) 

dQ T _ dg, 

dt dt dt 


total internal energy- 

specific enthalpy of 
fluid at the exit orifice 

total power input 


dM 

dt 


rate of fluid removal 
(negative for outflow) 


( 5 . 28 ) 

(5.29) 


(5.30) 


If a strictly conservative numerical method were used to integrate the exact 
conservation relations, equation (5.28) would he satisfied. A procedure for 
imposing equation (5.28) is developed in a later section. First, however, it 
is worthwhile to consider the corresponding energy conservation relation that 
applies to the approximate energy equation derived in chapter 2. 


Integration over the tank volume, application of Gauss' theorem to the 
energy equation derived in chapter 2, and use of M = pVip lead to the 
expression 


Me 


dT 
p dt 


*q t 

~~dt 


(K 



dM 

dt 


If a strictly conservative numerical method were used to integrate the equa- 
tions of chapter 2, the above equation would be satisfied. Conversely, overall 
energy conservation will result if this relation is imposed at the end of each 
time step by a uniform readjustment of the temperatures at the computa- 

tional grid points. A uniform readjustment is chosen because we have no infor- 
mation on the spatial location of the error. For this purpose h Q can be 
taken equal to c T and 

P „ dT _ /_ _ v 

Mc p dt dt (5.31) 


Setting h Q equal to CpT has the effect of removing the dM/dt term from 
equation (5*3l). The meaning of this choice is that we take the enthalpy at 
the exit orifice to be equal to the average specific enthalpy in the tank 
rather than attempting to specify an exact location of the exit orifice. Impos- 
ition of this relation at the end of each time step ensures conservation of 
energy and removes pressure excursions resulting from small but cumulative 
spurious internal energy sources arising from the non conservative numerical 
method described in chapter 3- A dimensionless temperature H is defined in 
equation (2.95a) and H ^ is evaluated at the computational grid points in the 
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numerical integration procedure of chapter 3 (eq. (3.2d)). The temperatures 
T^k at the computational grid points are related to the dimensionless tempera- 
tures by the relation 



+ T H 
air jk 


(5.32) 


where T R is a reference temperature that can be time dependent and T^^ is 
an arbitrary constant temperature difference. The uniform readjustment of the 
Tjk at each time step can be effected by a change in the reference temperature 
T r rather than modifying all elements of the computational matrix H-j^. For 
that purpose the volume average of equation (5*32) 


T = T R + T dif H (5.33) 

is needed. At the end of each time step the volume average H is computed to 
determine the necessary correction of Tp>. Integration of equation ( 5 • 31 ) over 
one time step yields 

- - 1 dQ T 

T ^ - T = tt ~ TT At 

n+l n Me dt 
P 

Substitution of equation (5.33) leads to the relation 


1 d0 T 


<*B> n+ , - <%>„ + WU - H n ) = Mc at 
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At 


from which it follows that the required correction = (Tr) 


n+l 


(TR> n is 


*] dQrp _ _ 

ATp = 7^ tv" At - T^.jH - H ) 

h Me dt dif n+l n 


(5.34) 


Correction to Account for Mean-Density Variations 

Additional nonphysical cumulative excursions of pressure can result from 
the use of approximate thermodynamics in the integration procedure for deter- 
mination of the motion and temperature of the fluid. If an exact compressible 
conservative method were used, equation ( 5 . 28 ) would be satisfied. Conversely, 
imposition of equation ( 5 . 28 ) by readjustment of dQp/dt will result in an 
overall energy balance and will restore the correct pressure variation with 
time. If this were done at every time step, the correction for spurious inter- 
nal sources (eq. (5-3^)) would be superfluous since it would be overridden. 
However, the computing time required for evaluating the necessary thermodynamic 
functions in their present form (described in ch. h) is not negligible. Since 
the correction to account for compressibility effects is relatively small for 
short time periods , computing time can be saved by imposing it only at every 
tenth time step. In that case, since the computing time required for imposing 
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the correction for spurious internal sources at every time step is negligible * 
both corrections are worthwhile. 


Substitution of equation (5-29) into (5*28) and integration over one time 
step yield 


( 55 ) 


n+i 


<p5>. 




— + h M 

at o dt , 


At 


This relation may not be satisfied if a nonconservative numerical method is 
used. It will be satisfied if a correction power dQ c /dt is added to the 
total input power dQrp/dt such that 




n+1 


* 0 ™ 


- h, 


dM 
o dt 


(5.35) 


It is convenient to add this small correction to the wall heat leak power 
dQ-^/dt discussed after equation (5.25). Thus heat is added (or subtracted) at 
the tank walls to allow for the difference between accurate thermodynamics and 
the approximate thermodynamics used in the integration procedure. Physically 
realistic pressure variations will result for cases in which the convection and 
conduction processes within the tank are adequately approximated by the two- 
dimensional integration procedure for a square tank described in chapters 2 and 
3 . 


Method for Including Effects of Tank Stretch and 
Variable Transport Properties 

Since the Apollo oxygen tanks are thin-walled pressure vessels, the tank 
volume depends on the pressure according to the relation 


dV„ 


V = V 
T T 


1 + 


af (p - p o> 


(5.36) 


where V T and p^ are the initial volume and pressure, respectively. Accord- 
ing to relations given in unpublished work of C. K. Forester, dV^/dp can be 
considered constant and is given by 


where 



3r(l - ct)V t 

o 

2b iE 

y 


(5.37) 


r/bi ratio of spherical .tank radius to wall thickness 
a Poisson’s ratio for tank wall material 
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initial tank volume 


v T 

o 

E Young’s modulus for tank wall material 

y 

The changes in tank volume due to changes in pressure are small, but they can 
have a large effect on the rate of pressure rise due to heating and on the 
potential pressure decay. To account for such effects it is necessary to 
replace equation (5.28) with 


dE 

dt 


dQ, 


— + h M _ 
dt o dt 


dV, 


_T 

dt 


which is a form of the first law of thermodynamics. This relation is applic- 
able to the entire volume of stratified fluid because pressure gradients are 
negligible. Differentiation of equation (5.36) with Vij , dV^/dp, and p held 
constant, and substitution in the last relation yield 0 u 

dE f^T dM 1 dV T dp 2 

dt dt o dt " 2 dp dt 

Integration of this equation leads to the expression 


E = E 


♦ /* ft -o !i) 4t - ^ (p2 - 


for the total internal energy in the tank. 


(5.38) 


To ensure conservation of energy in each time step a procedure similar to 
that leading to equation (5.2l) can he applied to the relation 


VA+l = (p N e N ) n + Ap rff 
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where is evaluated at n+1. Since pe = E/Vrp, substitution of equations 

(5.36) and (5.38) into the last relation yields 
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Similarly, use of p = M/Vrp and substitution of equation (5-36) into equation 
(5.21) yield 
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In this derivation F™ is evaluated at the (n+l)st integration step. If 
equations (5*39) and (5-40) are satisfied at each time step, overall conserva- 
tion of energy and mass is assured. At the end of an integration step all 
quantities appearing in these equations are known except Ap and AT R , the 
change in reference temperature imposed by the effect of tank stretch on the 
total energy. Note that equations (5*39) and ( 5.1+0) can be solved for Ap and 
AT r at each time step so that energy and mass are conserved. 


The temperatures Tj k at the computational grid points are related to the 
dimensionless temperatures computed in the integration procedure of chap- 

ter 3 by equation (5.32). Substitution of that equation into the expression 
for the distribution function F^ , equation (5.28) of chapter h yields 


b jtk 
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The Wj k are constant weighting values discussed previously. The T^ are an 
array of fixed temperatures with uniform spacing AT and T^^ is an arbitrary 
temperature difference. At the end of an integration step according to the 
method of chapter 3, the dimensionless temperatures are known. There- 

fore, all quantities needed to compute F^ for use in equations (5-39) and 
(5*1+0) are known except the reference temperature Tp. However, T R is known 
at the previous time step. Therefore, F^ is evaluated by means of the 
relation 
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( 5 . 2 + 1 ) 


F S “ VWJ + 4 % 


fV (T K ) n + 4T] ' 

AT 


where AT is the spacing between the T^, that is, AT = T^ +1 - T^. Substitut- 
ing equation (5-4l) into (5-39) and (5-40) yields a pair of equations too 
lengthy to be written out conveniently. However, the important point is that 
the resulting equations can be solved (by iteration) for Ap and AT^ at each 
time step. In this procedure no additional corrections for spurious internal 
sources or changes in mean density are required since energy and mass are 
already conserved. The structure of the relations in this section is such that 
the computed values of p and pe are automatically driven toward the correct 
values and significant drifts do not occur. 


At the end of each time step, and also at the start of the calculation, 
when p and are known a number of quantities needed for calculation at the 
next time step are evaluated at the fixed temperatures T^ and pressure p. 
Each element of the fixed double arrays 


’mn ~ p ^ p m ,t n^ 

(5.2+2) 

: MN = U( ' P M ,T N^ 

(5.2+3) 


used in these evaluations is computed and stored during the calculation the 
first time it is needed, using the thermodynamic relations of chapter 4. The 
stored values are used, if they have been previously computed, to minimize the 
computing time. It has been determined that a uniform spacing of 10 psi 
between the pjyj and 1° K spacing between the T^ provide adequate accuracy. 


The densities pjj 
evaluated according to 


and derivatives 
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where M is chosen such that p^ < p < p^ +1 . Additional derivatives needed in 
the calculation are 


_ P M,N+1 P MN 
" T lf + 1 - T W 


(5-2+7) 




116 



(5.U8) 
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Again M is chosen such that < p < p^ +1 and these derivatives are evalu- 

ated at all N so that the quantities 


[Qp/3T) p ] N 
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[(9e/9T) p ] N + 6 N p 


( 5 .^ 9 ) 

(5-50) 


can be computed and stored for all N. Arrays (kinematic viscosity) and 

kjj (thermal conductivity) are also computed and stored at the end of each time 
step for use in the next integration step. The relations used are 


V N = V(P N’ V 
k N = k ^ P N’ T N^ 


(5.51) 

(5.52) 


where the functions v(p,T) and k(p,T) are evaluated according to an unpub- 
lished recommendation of H. M. Roder. Equations ( 5 .U 9 ) through (5-52) indicate 
the evaluations used to impose variable coefficients in the equations of chap- 
ter 2. The value of N used at each computational grid point is chosen such 


that T w < T m < T. 


N 


Jk 


N+l ’ 


The tank stretch effect also modifies the relations for potential pressure 
decay listed after equation (5*9). The total mass of the collapsed state is 
unaffected, but the total energy and volume are given by 
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are also needed. Iterative solution of these relations is required to find 
the collapse pressure p 

At the beginning of an integration for the motion and temperature distri- 
bution of the fluid the initial pressure p Q , tank quantity, and dimensionless 
temperature distribution are specified. The tank quantity is defined to 

be the ratio of fluid mass in the tank to the fluid mass when the tank is full. 
Thus , the mass of oxygen in the tank is computed according to 

M = M full x Quantity (5.59) 

Then 


p = M/V 


(5.60) 


The distribution function F-^ is a function of T R ; see the discussion follow- 
ing equation (5-^0). Therefore, the relations 
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can be solved for T^ by iteration. Once Tp is known 
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E = PeV. 


(5.64) 


can be evaluated. 


Suppression of Nonphysical Excursions of Temperature and Vorticity 

At an early stage in this investigation it was noted that the numerical 
integration method described in chapter 3 results in small excursions of tem- 
perature outside the range of values imposed in the boundary conditions. This 
probably occurs because of the use of a coarse grid that does not adequately 
resolve the heater boundary layer. We have described methods for ensuring 
overall conservation of mass and energy. An additional remedy has been adopted 
to suppress nonphysical excursions of temperature. 

In simulating a heater, heater temperatures are computed at the beginning 
of an integration for each time step. Similarly, wall temperatures are 


.a 
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computed to account for. heater radiation absorbed by the walls and heat leak 
through the. walls. After these computations the highest and lowest tempera- 
tures in the entire temperature matrix are found. These limits are then 
imposed on the entire temperature matrix £,t the end of each integration step 
since they are valid physical limits. 

A similar type of nonphysical behavior can occur in the computation of 
vorticity due to the failure of the coarse grid to resolve boundary layers. 
Therefore, upper and lower limits were imposed on the vorticity matrix at the 
end of each time step in the same manner as for the temperature. When this was 
not done, much larger values of vorticity occurred at the. boundaries and pro- 
duced oscillations in the flow unless the time step was decreased. Imposition 
of vorticity limits stabilized the flow and allowed long computations to be 
made. Fixed limits of plus or minus four times the magnitude of the uniform 
vorticity corresponding to a rotation reversal starting at three revolutions 
per hour were imposed. These limits were chosen because they are beyond the 
extreme values observed at interior points in the flow for the cases investi- 
gated. There is no physical justification for vorticity limits as there is for 
temperature limits. There is a physical interpretation of the effect, however. 
Disturbances that extend over several computational grid points can be computed 
correctly, but small-scale disturbances that cannot be computed correctly with 
a coarse grid in any case are filtered out. The flow in thin boundary layers 
and near the corners of the rectangular enclosure will not be computed cor- 
rectly. However, the mixing action that occurs in the main body of the fluid 
for times up to several hours can be computed and should provide a conservative 
estimate of the mixing that would take place in a spherical tank without 
corners . 


Effect of Changes in Parameters 

Chapter 6 presents results from a series of calculations based on the 
methods of this and the preceding chapters. A IT x 17 computational grid was 



Figure 5.3— Effect of changes in parameters (in all 
cases, flow rate = 1 .5 lbm/hr, heat leak = 10 W, heater 
power = 1 15 W, rotation rate =0.4 rph except where 
noted, quantity = 90 percent except where noted). 


used in those calculations as well 
as fixed values of several other 
parameters. It is of interest to 
consider the effect of changing the 
total number of mesh points and 
other parameters. 

Figure 5.3 contains plots of 
potential pressure decay versus time 
during heater cycling. In all of 
the calculations represented the 
heater is turned on when the pres- 
sure falls below 870 psia due to gas 
removal and the heater is turned off 
when the pressure rises above 910 
psia. A flow rate of 1.5 lbm/hr, a 
heat leak rate of 10 watts , . and a 
heater power of 115 watts were 
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imposed in all of the calculations. A steady spacecraft rotation rate of 0.4 
revolutions per hour was used, except where noted, and the tank quantity 
(required in eq. (5-59)) was 90 percent except for two cases at 80 and 70 per- 
cent. The quantity decreases by about 1 percent in 2 hr when the mass-flow 
rate is 1.5 lbm/hr. 

The potential pressure decay rises when the heater is on and subsides when 
the heater is off, but there is a cumulative increase as the level of tempera- 
ture stratification increases. It can be seen in figure 5.3 that the heater is 
turned on and off a number of times in each of the calculations shown. The 
positions of the fourth, seventh, and tenth peaks are labeled in the figure as 
an indication of the effect of changes in parameters on the rate of pressure 
rise. Since the pressure fluctuates linearly between 870 and 910 nsia in all 
cases, the time required for 10 peaks in the potential pressure decay is 
inversely proportional to the rate of pressure rise. This is illustrated in 
sketch (a), which shows pressure and potential pressure decay versus time for 

two heating cycles. Note that the rate 
of pressure rise due to heating is 
reduced by the tank stretch effect. 

The amount of the reduction is 
inversely proportional to the ratio of 
times required to complete a given num- 
ber of cycles. Calculations with and 
without the tank stretch effect are 
compared in the plot at the top of fig- 
ure 5.3. The ratio of times required 
to reach the tenth peak indicates that 
the tank stretch effect decreases the 
rate of pressure rise by a factor of 
1.225 relative to the rate of rise for 
a rigid container. In unpublished work 
of the Propulsion and Power Division, 
NASA Manned Spacecraft Center, results 
were obtained showing that with uniform heating at 90 percent quantity the cor- 
responding ratio of rise rates is 1.563. Thus, we find the tank stretch cor- 
rection to be smaller for nonuniform stratified heating than it is for uniform 
heating. It can also be seen in the plots at the top of figure 5-3 that the 
tank stretch effect tends to suppress the cumulative rise in potential pressure 
decay but does not counteract it completely. Additional information on this 
point is supplied by comparison of potential pressure decays with and without 
tank stretch for an assumed linear temperature distribution across the tank. 

In a representative calculation at 90 percent quantity with a temperature 
spread of 20° K, the values of potential pressure decay computed were 135 psi 
with tank stretch and 208 psi for a rigid container. Thus, tank stretch 
becomes more effective at higher levels of stratification in counteracting the 
buildup of potential pressure decay that would be present in a rigid container. 

Effects of changes in other parameters can be noted in figure 5.3. Rota- 
tion rate refers to the steady rate of rotation of the spacecraft used in the 
calculation. In each case, the calculation was started with a uniform tempera- 
ture and with the fluid at rest relative to the rotating tank. In all calcula- 
tions the heater was located at an off-center position described more 



Sketch (a) 
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completely in chapter 6 (fig. 6.10). Five computational grid points forming 
a + were used, as heater elements . At the end of each time step in the- integra- 
tion, the stream function at these five points was set equal to the average of 
the values computed at the five points by the numerical method of chapter 3. 
This was done for the purpose of imposing zero velocity at the center of the 
heater. "Nonzero velocity at heater" in figure 5.3 refers to a calculation in 
which the stream function was not readjusted at heater elements. "Large heater 
volume" indicates a calculation in which 13 points were used as heater elements 
rather than five points. It can be seen in figure 5*3 that use of 13 points 
decreased the computed potential pressure decay by only a factor of about two. 
Most of the calculations in this report were computed using 17 * 17 uniformly 
spaced matrix points. In figure 5.3 results are given for a 33 x 33 matrix 
that are to be compared with the standard calculation represented by the solid 
curve in the top plot. For tank quantities less than 90 percent the rate of 
pressure rise and the cumulative rise in potential pressure decay are smaller 
as shown in the plots at the bottom of figure 5.3. 


REFERENCES 


1. Hirs chf elder , J.; Curtiss, C. F. ; and Bird, R. B.: Molecular Theory of 

Gases and Liquids. John Wiley & Sons, Inc., 195U, pp. 237s 250-252, 256. 

2. Borucki , W. J. : Absorption of Infrared Radiation by High-Density Oxygen. 

Letter to the Editor, J. Spacecraft and Rockets. 

3. Forester, C. K. : Numerical Integration of the Conservation Equations for 

the Storage and Expulsion of Single-Phase Cryogens . M.S. Thesis, Dept, 
of Mech. Engr., Univ. of Washington, 1969. 


121 




6. RESULTS FROM NUMERICAL COMPUTATIONS SIMULATING 


FLOWS IN AN APOLLO OXYGEN TANK 


Barrett S. Baldwin, Walter A. Reinhardt, and Yvonne S. Sheaf fer 


SUMMARY 


The effectiveness of vehicle maneuvering as a means for producing neces- 
sary mixing in the Apollo oxygen tanks is investigated. Results are presented 
from numerical simulation of flow in the tanks showing the effect of reversal 
of rotation after a prolonged period of rotation at a constant rate. Results 
are given from calculations corresponding to spin up of the vehicle after a 
prolonged period at zero rotation rate and near-zero acceleration. Both maneu- 
vers lead to a reduction in potential pressure decay by a factor of 2 or more. 
Photographs of a cathode ray display tube are presented illustrating the con- 
vection currents and cumulative buildup of potential pressure decay that result 
from intermittent operation of the heater. 


INTRODUCTION 


In this chapter, results are presented from a series of calculations uti- 
lizing the methods developed in preceding chapters. Preliminary results based 
on the van der Waals equations of state were presented at meetings on the 
Apollo oxygen system held at the Manned Spacecraft Center prior to the Apollo 
l4 flight. 

Following the failure of an oxygen tank during the Apollo 13 flight and 
the subsequent diagnosis of the cause, intensive efforts in the areas of 
design, development, analysis, and testing were conducted to assure the timely 
redesign of the Apollo cryogenic oxygen storage and supply system. After the 
Apollo l4 flight a symposium was organized at the NASA Manned Spacecraft Center 
for the purpose of disseminating the information acquired. Seventeen papers 
were presented under five categories: hardware development, stratification 

analyses, system models, test programs, and flight performance. The papers 
were collected in an unpublished document, MSC Cryogenics Symposium Papers. 

The relationship of the present investigation to the overall effort is outlined 
in chapter 1. 

The recommendation of the Apollo 13 Review Board that the mixing fans be 
removed from the oxygen tanks was adopted. This decision and the subsequent 
redesign studies were aided by data acquired from the Apollo 12 flight during 
which the fans were not used for an appreciable part of the flight. Analysis 
of that data was reported in unpublished work of C. K. Forester, D. D. Rule, 
and H. W. Patterson. Of particular interest was a period from 4 to 8 hr after 
the launch during which the fans were not operated. During most of this time 
very little spacecraft maneuvering occurred, and the effective acceleration 
field was of order 10 “ 7 earth gravity. It was found that a number of pressure 
decays took place that were correlated with accelerations due to maneuvering, 
which caused spikes of acceleration of order 10“ 4 g lasting for several 
seconds. After a prolonged period of low acceleration a pressure decay of l40 
psi was observed at the time of an acceleration spike that reached a level of 
10 “ 3 g. A stratification analysis based on a two-dimensional square-tank simu- 
lation reproduced these effects reasonably well. 
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The analyses of Apollo 12 data indicated that no serious stratification 
problems were to be expected from operation of the oxygen tanks with the mixing 
fans removed. It was concluded that normal vehicle maneuvering could be relied 
on to produce the necessary mixing. Nevertheless, several investigations, were 
made of emergency procedures that could be used if the performance fell below 
an acceptable level. The question was raised whether reversal of the rotation 
rate starting at three revolutions per hour (rph) would provide useful mixing 
action in the oxygen tanks. An analysis of the stratification using a two- 
dimensional model for this problem is difficult because prolonged motion of the 
fluid leads to the requirement of a small time step to achieve stability in the 
numerical method. Therefore, approximate methods were used initially to reduce 
the required computing time. From our preliminary results based on the van der 
Waals equations it was concluded that spacecraft rotation reversals would 
indeed provide effective backup means for mixing the oxygen if normally occur- 
ring maneuvers proved to be insufficient. That conclusion is of continuing 
interest for future Apollo flights in which vehicle maneuvering is relied on to 
produce the necessary mixing in the oxygen tanks. 

In this chapter we present results based on the accurate thermodynamic 
relations developed in chapter 4. The effects of tank stretch analyzed in 
chapter 5 are accounted for. The equations for the fluid motion and convection 
of heat used in the analysis are derived in chapter 2. The numerical method 
employed in the integration of the fluid mechanical equations is described in 
chapter 3. Special methods for suppressing cumulative errors that were used in 
the calculations are described in chapter 5* Our preliminary conclusions on 
the mixing effectiveness of rotational spacecraft maneuvers are verified and 
put on a firmer basis by the results in this chapter. The effects of changes 
in the rate of spacecraft rotation on additional types of stratification are 
investigated. Calculations showing the buildup of potential pressure decay 
(the pressure decay that would result from complete adiabatic mixing) during 
heater cycling are presented. 


MIXING EFFECTIVENESS OF ROTATION REVERSAL 


Calculations were made in which the vehicle was taken to be rotating ini- 
tially at 3 rph. A series of initial stratified states were imposed with the 
temperature varying linearly across the tank. The hot fluid was placed in the 
stable position toward the center of rotation, which was outside the tank. 

Such stratification can be expected to develop after many heater cycles while 
operating with a steady vehicle rotation rate. In the absence of other vehicle 
maneuvering and with the heater turned off, the calculations show a very slow 
decrease in potential pressure decay and no motion of the fluid relative to the 
tank. The decrease in potential pressure decay in this case is due to conduc- 
tion arising from the mild temperature gradient. When the direction of rota- 
tion is abruptly reversed, however, a swirling motion of the fluid ensues, 
leading to mixing and enhanced temperature gradients. Figure 6.1 depicts the 
velocity field in the flow that results. It comprises photographs of a cathode 
ray display tube on which were plotted the velocity vectors at the computa- 
tional grid points. Sketch (a) shows the dimensions and orientation of the 
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Figure 6.1 - Velocity vectors after rotation reversal. Sketch (a) 


two-dimensional square-tank model. As indicated, the tank rotates counter- 
clockwise in the plane of the page about a center located at the center of the 
spacecraft. However, in the calculations a coordinate system is used that is 
fixed with respect to the tank and the center of rotation remains fixed at a 
position 1.5 tank diameters below the center of the tank. In all the photo- 
graphs of figure 6.1 and later figures, the center of rotation is thus below 
the tank as indicated in the sketch. The arrow on the right in the photographs 
represents the position of a fixed star, such as the sun, with respect to the 
rotating tank. Thus when the rotation rate is steady, the arrow provides an 
indication of the passage of time. This feature was useful in viewing the 
cathode ray tube during calculations or in motion pictures of the displays. 

The upper left photograph in figure 6.1 shows the fluid velocity vectors in the 
tank immediately after rotation reversal. The magnitude of the velocity near 
the tank boundaries is about 0.02 ft sec"* 1 . The clockwise swirling motion 
shown results in part from the rotational inertia of the fluid, which tends to 
retain the motion it possessed before the rotation reversal of the vehicle. A 
lateral acceleration, present during the reversal, acts differentially on the 
stratified layers resulting in a significant contribution to the swirling 
motion. The other photographs in figure 6.1 show the velocity vectors at later 
times. The off-center swirl develops as a result of the stratification present 
and moves continuously in the clockwise direction, which is also the direction 
of motion of the fluid. 

The graph at the bottom of each photograph in figure 6.1 is a plot of the 
potential pressure decay versus time from the beginning of the calculation. 

The scale is automatically decreased when the plot becomes overextended. In 
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the lower right photograph the time scale has heen decreased hy a factor of 2 
relative to that in the other photographs . 

Figure 6.2 shows the evolution of temperature distribution in the tank. 

For temperatures greater than the average, the deviations from the mean temper- 
ature are represented by vertical lines. As an aid to visualization, for tem- 
peratures less than average the deviations from the mean are shown as horizon- 
tal lines. In either case, the length of the lines indicates the relative 
magnitude of the temperature deviation from the mean at each computational grid 
point. The upper left photograph in figure 6.2 shows the initial assumed 
linear temperature distribution. The hot gas is toward the center of rotation. 
The remaining photographs show temperature distributions at later times after 
mixing has resulted from the swirling motion. 

In figure 6.3 thermal stratification is viewed differently. The tempera- 
ture distribution function described in chapter 4 (following eq. (4.24)) and in 
chapter 5 (preceding eq. (5*13)) is displayed in the form of histograms that 
indicate the number of computational grid points with temperatures within 1° K 
intervals. The unfilled histogram represents the initial linear temperature 
distribution. Gaps appear between the bars due to the discretization used in 
the numerical method. That is, every other 1° K interval did not happen to 
contain temperatures occurring at the computational grid points in the imposed 
linear temperature variation across the tank. The bars at 130° and 170° K are 
half the size of the others because the computational grid points at the tank 
boundaries are given half the weight of interior points, which represent twice 
as much fluid volume. During a calculation the temperatures change due to con- 
vection and conduction of heat. The heights of the bars change and the gaps 
are filled in. The shaded histogram in figure 6.3 shows the temperature distri- 
bution 40 min after a rotation reversal maneuver when considerable mixing has 
taken place and the steep temperature gradients that developed have reduced the 
temperature deviations from the average . If complete mixing were to take 
place, the temperatures at all computational grid points would fall within the 
same 1° K interval and the histogram would become a single bar at that tempera- 
ture. As a supplement to the potential pressure decay, histograms of this type 
provide additional quantitative information on the stratification present in 
the tank. 

In a preliminary unpublished paper we have shown that under certain condi- 
tions the decrease in potential pressure decay due to a rotation reversal 
according to a calculation based on the van der Waals equations is in rough 
agreement with results based on more exact thermodynamic relations. However, 
such agreement is not obtained in all cases of interest, as illustrated by a 
comparison of figures 6.4 and 6.5. Figure 6.4 shows plots of potential pres- 
sure decay versus the magnitude of linear temperature variations for several 
tank pressures according to the van der Waals equations. Figure 6.5 shows 
results for the same conditions based on the accurate thermodynamic relations 
of Stewart as described in chapter 4. It is easily seen that the van der Waals 
results disagree grossly with those based on Stewart* s equations at both low 
and high levels of stratification. Our preliminary conclusions on the mixing 
effectiveness of a spacecraft rotation reversal although qualitatively correct 
are put on a firmer basis by the calculations utilizing accurate thermodynamics 
in this chapter. 
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Figure 6.2.— Temperature distributions after rotation 
reversal. 
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Figure 6.3.— Temperature distribution functions. 
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temperature variations; van der Waals equations. 



Figure 6.5.— Potential pressure decay for linear 
temperature variations; Stewart’s equations. 


The sharp bends in the curves contained in figure 6.5 are of interest. 

The flattening out takes place when a level of stratification is reached for 
which the collapsed state that would result from complete adiabatic mixing con- 
tains a two-phase mixture of liquid and vapor. Many other interesting aspects 
of the behavior of cryogenic oxygen are illustrated in the thermodynamic prop- 
erty plots presented in chapter k. 


Initial pressure 900psio 
Mean temperature I50°K 


Initial rotation rate 3rph Courant No„= 0.4 



Figure 6.6 — Potential pressure decay after rotation reversal, 
equation (3.28) can be written 


Figure 6.6 shows plots of 
potential pressure decay versus 
time after rotation reversal 
for four initial levels of 
stratification. The starting 
values of potential pressure 
decay in these plots can be 
obtained from figure 6.5* 
Decreases in potential pressure 
decay by a factor of about 2 or 
more result from a change in 
rate of spacecraft rotation 
from 3 rph to -3 rph. 

The Courant number applic- 
able to the numerical integra- 
tion procedure used in the 
present calculations given in 


..M 


At 


Ax 

where |u| is the maximum speed present in the flow. At the time step, and 
Ax the spacing between computational grid points in both x and y directions. 
As discussed in chapter 3, the Courant number is a measure of the stability and 



accuracy of the numerical integration procedure. Most of the calculations in 
this report were made with a value of a = 0.8. Case I of figure 6.6 was also 
computed with a = 0.U (and smaller time steps according to eq. (6.1)) as indi- 
cated by the dotted line. Comparison of the dotted and solid line curves for 
case I in figure 6.6 shows that the predicted variation of potential pressure 
decay with time is insensitive to such a change in the Courant number. Appre- 
ciable variations in the flow variables extend over several computational grid 
points except in the boundary layer at the wall, which remains relatively thin 
during computation times considered here. Therefore, changes in the boundary 
layer do not significantly affect the evolution of the potential pressure 
decay, which depends on conditions in the main bulk of the fluid, in the rota- 
tion reversal problem. For this reason, the results are insensitive to a 
change in grid spacing from a 17 * 17 matrix to a 33 * 33 matrix. To save com- 
puting time, the latter computations were also made with a Courant number 
a = 0.8. 


MIXING EFFECTIVENESS OF 'SPINUP AFTER ATTITUDE HOLD 


In the previous section, initial stratified states were considered that 
could be expected to result after many heater cycles with the spacecraft rotat- 
ing at a steady rate (the so-called "passive thermal control mode"). Another 
case of interest is the type of stratification that would result after many 
heater cycles in a nonrotating state (attitude hold mode). In the absence of 
vehicle maneuvering, no convection currents would develop and the heat from the 
heater can spread into the gas only by conduction. A localized hot spot around 
the heater somewhat diffused by conduction and essentially zero fluid velocity 
are to be expected in this case. It is of interest to determine the mixing 
effectiveness of a spinup to a steady rotating state from such an initial strat- 
ified state. In the coordinate system fixed with respect to the tank, changes 
in rotation rate cause a rotating motion of the fluid. The velocity field that 
occurs as a result of spinup is similar to that shown in figure 6.1. figure 
6.7 illustrates the effect of spinup on the temperature distribution. It 
should be recalled that positive temperature deviations (from the mean tempera- 
ture) are shown as vertical line segments and negative deviations are shown as 
horizontal line segments. The upper left photograph shows the initial assumed 
distribution. Subsequent distortions and dissipation resulting from the 
swirling motion are shown in the remaining photographs . 

Figure 6.8 is a plot of potential pressure decay versus time after spinup. 
For comparison, a plot is included showing the very slow decrease in potential 
pressure decay that occurs as a result of conduction when the attitude hold 
condition is maintained. Again it is found that a change in vehicle rotation 
rate provides effective mixing action. 

HEATER CYCLE 


A pressure sensor is installed in the Apollo oxygen tanks 
automatically turn the heater on when the pressure drops below 


and is used to 
a lower limit of 
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Figure 6.7.— Temperature distributions after spinup. Figure 6.8 — Reduction in potential pressure decay due 

to spinup. 


about 870 psia as a result of gas removal or pressure decay due to mixing. The 
heater is automatically turned off when the pressure rises above an upper limit 
of about 930 psia due to heating. As discussed in chapter 5 our machine pro- 
gram contains options for simulating the heater in either of two ways . One 
method utilizes a segment of the tank wall as a heater. In the other method, 
selected internal computational grid points can be utilized as heater elements 
to more nearly correspond to the actual position of the heater in an Apollo 
oxygen tank. 

The photographs of figure 6.9 show the velocity and temperature distribu- 
tions produced by cyclic operation of the two types of heater. The plots at 
the bottom of each photograph show the variations in potential pressure decay 
that have taken place since the beginnings of the calculations. When the 
heater is on, the potential pressure decay rises, and when the heater is off, 
it subsides, but there is a cumulative increase. About 2 hr (real time) have 
elapsed since the start of the calculation at the time of the photographs . 
During that time about six heater cycles have been completed. Additional 
information on the parameters of the calculations is given in the discussion of 
subsequent figures. 

The upper left photograph in figure 6.9 shows the velocity vectors result- 
ing from operation of a heater on the left wall. The upper right photograph 
shows the temperatures at the computational grid points. The temperature of 
hot fluid is represented by vertical line segments and the temperature of fluid 
cooler than average is indicated by horizontal lines. In either case, the 
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Figure 6.9.— Velocity and temperature distributions 
resulting from heater operation (0.4 rph). 


heater (0.4 rph). 

length of the line segment repre- 
sents the deviation of the temper- 
ature from a reference value. 

Each time the heater is turned 
off, the reference temperature 
used in the display is readjusted 
to the value of the mean tempera- 
ture in the tank. This is done so 
that hot fluid can he easily dis- 
tinguished from cold fluid in the 
displays shown. For comparison , 
the velocity and temperature dis- 
tributions arising from operation 
of an internal heater are shown in 
the lower two photographs of 
figure 6-9- 

Figures 6.10 and 6.11 contain 
plots of pressure and potential 
pressure decay for the two heater 
positions. A steady spacecraft 
rotation rate of 0.U rph was 
imposed. Additional parameters of 
the calculations are indicated on 
the figures. The notation 

QUANTITY = 90 percent indicates that the tank contains 90 percent of the amount 
of fluid it contained when full. FLOW RATE =1.5 lbm/hr refers to the rate of 
fluid removal from the tank. HEAT LEAK - 10 W indicates that heat is leaking 
into the fluid from the exterior of the tank at the rate of 10 W. When the 
heater is on it produces heat at the rate of 115 W. The pressure rises due to 
heating when the heater is on and falls due to gas removal when the heater is 
off. Upper and lower pressure limits were imposed corresponding to data from 



Figure 6.11.— Pressure cycles for square tank with 
internal heater (0.4 rph). 
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portions of the Apollo lb flight. The restilts in figure 6.10 were generated 
using a wall heater * and the calculation represented by figure 6.11 was identi- 
cal except for the use of an internal heater position. Comparison of the two 
figures shows that there are no gross differences arising from the position of 
the heater for the conditions of these calculations. 


Figure 6.12 contains photographs showing velocity and temperature fields 
produced by prolonged heater cycling of an internal heater with the spacecraft 
rotating at a steady rate of 0.1| rph. The upper photographs were taken at a 
time corresponding to about 8 hr after the start of the calculation and the 
time of the lower photograph is 10 hr. The velocities remained small (less 
than 10'” 4 ft /sec) so that the stability criteria of the numerical method 
allowed a large time step (0.5 min) in the integration. The total computing 
time on an IBM 360-67 computer was b6 min. The photographs of temperature dis- 
tribution on the right of figure 6.12 are somewhat similar to the assumed ini- 
tial linear distributions 1 used for previous calculations in this article. With 
prolonged heater cycling during steady spacecraft rotation, it is seen that 
there is a tendency for the hot gas (represented by vertical line segments) to 
collect at the bottom of the tank in the direction of the center of rotation. 

On the plots of potential pressure decay at the bottom of each photograph in 
figure 6.12, in the lower two frames the time scale has been halved relative to 
that in the upper two frames. These plots show that the cumulative rise in 
potential pressure decay does not continue indefinitely when the spacecraft is 
rotating. Instead a plateau is reached when the velocity past the heater has 
built up sufficiently that the incoming heat is convected to other regions of 
the tank rather than continuing to accumulate in a small region near the 
heater. As this natural convection effect continues to grow, there is a 

velocity Temperature decrease in the level of the 



Figure 6.12.— Velocity and temperature distributions 
after prolonged heater cycling at 0.4 rph. 


potential pressure decay. In the 
lower two frames of figure 6.12, 
it can be seen that near the end 
of the calculation the cumulative 
increase in potential pressure 
decay has resumed. This is prob- 
ably due to the accumulation of 
hot gas at the bottom, which 
tends to counteract the driving 
force of the natural convection 
and reduce the velocity of flow 
past the heater. The highest 
level of potential pressure decay 
reached in this calculation was 
about 80 psi. 

Figure 6.13 shows the effect 
of a change in spacecraft rota- 
tion rate on the cumulative rise 
in potential pressure decay. In 
this calculation the spacecraft 
was not rotating initially. The 
cyclic heater operation resulted 
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in a rise of potential pressure decay 
to a level above 100 psi. At a real 
time 260 min after the start of the 
calculation, a spinup to 3 rph was ini- 
tiated and completed in about 10 min. 
The rotation rate was held steady 
thereafter. Figure 6.13 shows that the 
change in spacecraft rotation rate 
resulted in an appreciable decrease in 
the potential pressure decay. The 
forced convection introduced by the 
0 40 80 i 2 o 160 200 240 280 320 maneuver entails larger fluid veloci- 

T ime, mm 

ties than those arising from natural 
Figure 6.13. — Effect of spinup on cumulative rise in convection in the previous case con- 
potential pressure decay. sidered. Furthermore , these larger 

velocities persist rather than being 
suppressed by accumulations of hot 
fluid on the side of the tank toward 
the center of rotation. Consequently, the potential pressure decay would be 
expected to continue to decrease if the calculation had been continued to 
longer times. The machine computing time on this run was 64 min. About two- 
thirds of that time was spent in the part of the calculation after the change 
in rotation rate. As a result of the maneuver the velocities became relatively 
large (~0.005 ft/sec), and the size of time step required for numerical stabil- 
ity decreased by more than an order of magnitude. 

DISCUSSION 


It is appropriate to discuss the probable effect of several approximations 
that were employed in this investigation. Two-dimensional simulation of the 
flow in a spherical tank was necessitated by the limited speed and storage 
capacity of the currently available computational facilities . If a highly 
stratified state is present in the tank with an appreciable fraction of the 
oxygen at a temperature 10° or more above the average temperature, a small tem- 
perature gradient will develop. Due to the low thermal diffusivity of cryo- 
genic oxygen, such a nonuniform condition can persist for many hours if the 
flow velocity is negligible. The fluid velocity induced by spacecraft maneu- 
vering will lead to penetration of the hot fluid by cold streams and lengthen 
the boundary between hot and cold regions. The enhanced temperature gradients 
and lengthened regions of heat exchange promote a more rapid trend toward a 
uniform temperature. The extent to which this occurs in a two-dimensional sim. r 
ulation should be indicative of the same trend in the actual three-dimensional 
flow if the driving force is in the direction of the two dimensions considered. 
There is no apparent reason to expect a large order-of-magnitude error in the 
estimate of mixing effectiveness from the two-dimensional calculations. 

It was not possible to adequately resolve the thermal and viscous boundary 
layers that would occur in a real flow. Estimates based on the viscous flow 
over a flat plate indicate that the boundary layer at the tank wall due to 
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changes in rate of spacecraft rotation would be laminar. Also the boundary 
layer would not encompass an appreciable fraction of the fluid in the tank for 
periods up to 2 hr after a change in rate of rotation, if the fluid were ini- 
tially at rest relative to the tank. During such times, the calculations show- 
appreciable mixing in the main body of the fluid due to interaction of the 
acceleration field with the nonuniform density distribution in the stratified 
field. Therefore, it is not expected that accurate treatment of the boundary 
layer at the wall would appreciably alter our predictions of the response of 
the potential pressure decay to spacecraft maneuvers. The thermal and viscous 
boundary layer that develops at the heater was also not adequately resolved in 
our calculations. However, the nonuniform temperature distributions produced 
by the approximate heater simulation provide worthwhile tests of the mixing 
effectiveness of spacecraft maneuvers . 

An objection can perhaps be raised to the use of a coarse computational 
grid spacing. Physical instabilities may develop in the boundary between hot 
and cold regions of a moving fluid. Generally, it is not possible to distin- 
guish such an occurrence from a numerical instability. Special methods that 
have been used to control the numerical stability may therefore have eliminated 
physically real effects. In that case, the calculations should provide a con- 
servative estimate of the mixing effectiveness of a spacecraft maneuver. 

Because of the previous existence of numerical methods for computation of 
fluid flows in rectangular enclosures, it was expedient to consider first the 
flow in a rotating tank of square cross section. The coarse grid spacing and 
the special methods employed to control stability may have obscured physically 
real effects that would occur near the corners in an actual fluid flow in a 
rotating enclosure of square cross section. The error from such effects would 
be expected to grow with time and could be appreciable for long computation 
times. For the present purpose, omission of such corner effects is not objec- 
tionable since our purpose is to simulate the flow in a spherical tank without 
corners. We have developed a machine program for computing the two-dimensional 
flow in a rotating tank of circular cross section. Preliminary results from 
that program were presented in unpublished reports at the MSC Cryogenics Sym- 
posium. No gross differences were found in predictions of the mixing effec- 
tiveness of changes in spacecraft rotation rate based on calculations for 
square or circular cylindrical tank geometries. 

The present investigation has been confined to large tank quantities 
(large values of the ratio of fluid mass to the mass when the tank is full) . 
Unpublished work of C. K. Forester, D. D. Rule, and H. W. Patterson has shown 
that the levels of potential pressure decay to be expected decrease rapidly 
with decreasing tank quantity for values of tank quantity below about 70 per- 
cent. This occurs because of the departure of the fluid state from the prox- 
imity of the critical point as fluid is removed from the tank and the fluid 
density decreases. At a tank quantity of 10 percent the properties of oxygen 
are approximately those of a thermally and calorically perfect gas for which 
the potential pressure decay is zero regardless of the degree of stratifica- 
tion. For this reason, from the operational point of view, the response of the 
fluid states to spacecraft maneuvers is of greater interest when the tank 
quantity is large. 



CONCLUDING REMARKS 


The effectiveness of rotational spacecraft maneuvers as a means for pro- 
ducing necessary mixing in the Apollo oxygen tanks has been investigated. 
Photographs of a Cathode ray display tube have been presented illustrating the 
convection currents resulting from changes in rate of rotation, and from natural 
convection. Results have been presented showing the cumulative buildup of 
potential pressure decay that occurs due to intermittent operation of the 
heater used to maintain the pressure in the tank during removal of fluid. 

There are several significant contributions from the present investigation 
toward a better understanding of conditions in the Apollo oxygen tanks during 
flight. We find that the mixing due to interaction of fluid density gradients 
with the acceleration field from rotational spacecraft maneuvers takes place in 
times that are short compared to those required for mixing due to viscous 
effects. For example, unpublished results from water tank simulation experi- 
ments of J. F. Lands, Jr., and R. C. Ried, Jr., indicate that times of order 4 
to 20 hr are required for the effects of a spacecraft rotation reversal to 
spread throughout the bulk of the fluid. In contrast, our results, which 
include the effects of density gradients, show considerable mixing throughout 
the fluid volume in less than 1 hr. For reasons previously discussed, our cal- 
culations represent a conservative estimate of the amount of mixing to be 
expected in actual flight. 

The levels of potential pressure decay to be anticipated according to our 
calculations are in reasonable agreement with previous estimates from Apollo 12 
data and stratification analyses. We concur in the conclusion that dangerous 
levels of stratification will not occur with mass flow rates and acceleration 
fields corresponding to past and anticipated Apollo operational practices. If 
as a result of departures from usual practice, an unacceptable level of strati- 
fication is believed to exist during flight, we find that return to a safe 
level can be achieved by mild changes in the rate of spacecraft rotation. A 
change in rotation rate of three revolutions per hour produces a reduction in 
potential pressure decay by a factor of 2 or more in 1 hr. Again, this repre- 
sents a conservative estimate. 

In the event of future space-flight applications that require scaling of 
the Apollo design with respect to time or size, the computational methods of 
this report are applicable for estimating the levels of stratification and the 
response to spacecraft maneuvers to be anticipated. 


'' NASA-Langley, 1972 12 A- 4224 
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